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Abstract 

This paper presents an overview and introduction to Smoothed Particle Hydrodynamics and Magnetohydro- 
dynamics in theory and in practice. Firstly, we give a basic grounding in the fundamentals of SPH, showing 
how the equations of motion and energy can be self-consistently derived from the density estimate. We 
then show how to interpret these equations using the basic SPH interpolation formulae and highlight the 
subtle difference in approach between SPH and other particle methods. In doing so, we also critique several 
'urban myths' regarding SPH, in particular the idea that one can simply increase the 'neighbour number' 
more slowly than the total number of particles in order to obtain convergence. We also discuss the origin of 
numerical instabilities such as the pairing and tensile instabilities. Finally, we give practical advice on how 
to resolve three of the main issues with SPMHD: removing the tensile instability, formulating dissipative 
terms for MHD shocks and enforcing the divergence constraint on the particles, and we give the current 
status of developments in this area. Accompanying the paper is the first public release of the ndspmhd SPH 
code, a 1, 2 and 3 dimensional code designed as a testbed for SPH/SPMHD algorithms that can be used to 
test many of the ideas and used to run all of the numerical examples contained in the paper. 

Keywords: particle methods; hydrodynamics; Smoothed Particle Hydrodynamics; Magnetohydrodynamics 
(MHD); astrophysics 



1. Introduction 



Smoothed Particle Hydrodynamics (SPH), originally formulated by Lucy (1977) and Gingold and Mon- 



aghan (19771, is by now very widely used for many diverse applications in astrophysics, geophysics, engi- 



neering and in the film and computer games industry. Whilst numerous excellent reviews already exist (e.g. 



Monaghan 1992 2005 Price 2004 Rosswog 2009 1 , there remain - particularly in the astrophysical domain 



- some widespread misconceptions about its use, and more importantly, its fundamental basis. 

Our aim in this - mainly pedagogical - review is therefore not to provide a comprehensive survey of 
SPH applications, nor the implementation of particular physical models, but to re-address the fundamentals 
about why and how the method works, and to give practical guidance on how to formulate general SPH 
algorithms and avoid some of the common pitfalls in using SPH. Since such an understanding is critical to the 
development of robust and accurate methods for Magnetohydrodynamics (MHD) in SPH (hereafter referred 
to as "Smoothed Particle Magnetohydrodynamics", SPMHD), this will lead us naturally on to review the 
background and current status in this area - particularly relevant given the importance of MHD in most, 
if not all, astrophysical problems. Whilst the paper is written with an astrophysical flavour in mind (given 



*This review and associated material germinated as lectures and tutorials given as part of the ASTROSIM summer school 
on computational astrophysics held during July 2010 in Torun, Poland. A video of the original lectures can be viewred online 
at http://supercomputiiig.astri.umk.pl/ The paper also presents an updated version of much of the otherwise unpublished 
material in my PhD thesis (Price 2004). 
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Figure 1: Computing a continuous density field from a collection of point mass particles, a) In particle-mesh methods (left 
panel) the density is computed by interpolating the mass to a grid (or simply dividing the mass by the volume). However, this 
tends to over/under resolve clustered/sparse regions, b) An alternative not requiring a mesh is to construct a local volume 
around the sampling point, solving the clustering problem by scaling the sample volume according to the local number density 
of particles, c) This panel shows the approach adopted in SPH, where the density is computed via a weighted sum over 
neighbouring particles, with the weight decreasing with distance from the sample point according to a scale factor h. 



the topical issue of JCP for which it is written), the principles are general and thus are applicable in any of 
the areas in which SPH is applied. 

Finally, alongside this article I have released a public version of my ndspmhd SPH/SPMHD code, along 
with a set of easy-to-follow numerical exercises - consisting of setup and input files for the code and step-by- 
step instructions for each problem in 1, 2 and 3 dimensions ~ the problems themselves having been chosen 
to illustrate many of the theoretical points in this paper. Indeed, ndspmhd has been used to compute all 
of the test problems and examples shown. The hope is that this will become a useful resourccj^ not only for 
advanced researchers but also for students embarking on an SPH-based research topic. 



2. The foundations of SPH: Calculating density 

The usual introductory lines on SPH refer to it as a "Lagrangian particle method for solving the equations 
of hydrodynamics". However, SPH starts with a basis much more fundamental than that, as the answer to 
the following question: 

How does one compute the density from an arbitrary distribution of point mass particles? 

This problem arises in many areas other than hydrodynamics, for example in obtaining the solution to 
Poisson's equation for the gravitational field = 47rG'p(r) when a (continuous) density field is represented 
by a collection of point masses. 



2.1. Approaches to computing the density 

Three common approaches are summarised in Fig. [l] Perhaps the most straightforward (Fig. [T^) is to 
construct a mesh of some sort and divide the mass in each cell by the volume. This basic approach forms the 
basis of hybrid particle-mesh methods such as Marker-In-Cell (e.g. Harlow and Welch 1965| ) and Particle-In- 
Cell (Hockney and Eastwood 1981 ) schemes, where one can further improve the density estimate using any 
of the standard particle-cell interpolation methods, such as Cloud-In-Cell (CIC), Triangular-Shaped-Cloud 
(TSC) etc. However there are clear limitations - firstly that a fixed mesh will inevitably over/under-sample 
dense/sparse regions (respectively) when the mass distribution is highly clusteretQ and secondly a loss of 
accuracy, speed and consistency because of the need to interpolate both to/from the particles, for example 
to compute forces. 



^NDSPMHD is available from http://users.inonash.edu.au/~dprice/SPH/ Note that we do not advocate the use of NDSPMHD 
as a "performance" code in 3D, since it is not designed for this purpose and excellent parallel 3D codes already exist (such as 
the GADGET code by |Springel|2005| . Rather it is meant as a testbed for algorithmic experimentation and understanding. 

^More recently, this problem has been addressed by the use of adaptively refined meshes to calculate the density field (e.g. 
|Couchman"||1991 1. 
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The second approach (Fig. [T|d) is to remove the mesh entirely and instead calculate the density based on 
a local sampling of the mass distribution, for example in a sphere centred on the location of the sampling 
point (which may or may not be the location of a particle itself). The most basic scheme would be to divide 
the total mass by the sampling volume, i.e., 

The problem of resolving clustered/sparse regions can be easily addressed in this method by adjusting the 
size of the sampling volume according to the local number density of sampling points, for example by 
computing with a fixed "number of neighbours" for each particle - as shown in the Figure. However, this 
leads to a very noisy estimate, since the density estimate will be very sensitive to whether a distant particle 
on the edge of the volume is "in" or "out" of the estimate (with Sp cx l/N^eigh for equal mass particles). This 
leads naturally to the idea that one should progressively down-weight the contributions from neighbouring 
particles as their relative distance increases, in order that changes in distant particles have a progressively 
smaller influence on the local estimate (that is, the density estimate is smoothed). 

2.2. The SPH density estimator 

This third approach forms the basis of SPH and is shown in Fig. [T|:: Here the density is computed using 
a weighted summation over nearby particles, given by 

p(r)= J2 mbW{r^n,h), (2) 

b=l 

where W is an (as yet unspecified) weight function with dimensions of inverse volume and h is a, scale pa- 
rameter determining the rate of fall-off of as a function of the particle spacing (also yet to be determined). 
Conservation of total mass J pdV — X^h^]"^' ™6 implies a normalisation condition on W given by 

/ Ty(r'~r,,/i)dV' = l. (3) 
Jv 

The accuracy of the density estimate then rests on the choice of a sufficiently good weight function 
(hereafter referred to as the smoothing kernel). Elementary considerations suggest that a good density 
kernel should have at least the following properties: 

1. A weighting that is positive, decreases monotonically with relative distance and has smooth derivatives; 

2. Symmetry with respect to (r — r') - i.e., W{r' — r,h) = W{\r' — r\,h); and 

3. A flat central portion so the density estimate is not strongly affected by a small change in position of 
a near neighbour. 

A natural choice that satisfies all of the above properties is the Gaussian: 



W{r-r\h) = ^cxp 



l\2 



(r-r') 



h 



2 



(4) 



where d refers to the number of spatial dimensions and cr is a normalisation factor given by cr = [1/ y^, I/tt, 1/ (tt^ 
in [1,2,3] dimensions. The Gaussian satisfies condition 1 particularly well since it is infinitely smooth (dif- 
ferentiable) - and gives in practice an excellent density estimate. However it has the practical disadvantage 
of requiring interaction with all of the particles in the domain [with computational cost of 0{N'^) if comput- 
ing the density at the particle locations], despite the fact that the relative contribution from neighbouring 
particles quickly becomes negligible with increasing distance. Thus in practice it is better to use a kernel 
that is Gaussian- like in shape but truncated at a finite radius (e.g. a few times the scale length, h). Using 
kernels with such "compact support" means a much more efficient density evaluation, since the cost scales 
like 0{NneighN), but inevitably leads to a more noisy density estimate since one is more sensitive to small 
changes in the local distribution. 
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Figure 2; The M4 (cubic, truncated at 2h), M5 (quartic, truncated at 2.5h) and Mg (quintic, truncated at 3fe) |Schoenberg| 
l |1946| B-spline kernel functions (solid lines) and their first (long-dashed) and second (short-dashed) derivatives, compared to 
the Gaussian (right panel and dotted lines in other panels). Notice that although the "number of neighbours" increases in 
the M5 and Me functions compared to the cubic spline, the smoothing scale h retains the same meaning with respect to the 
Gaussian. Thus, using the higher order B-splines is a way to increase the smoothness of the kernel summations without altering 
the resolution length, and is very different to simply increasing the number of neighbours under the cubic spline. 



2.3. Kernel functions with compact support 

There are many kernel functions which fit this bill. The most well-used (for SPH at least) are the 
Schoenberg (1946) B-spline functions (Monaghan and Lattanziol 1985 Monaghanl 1985 2005), generated 



as the Fourier transform 



M„(x,/i) 



1 r°° 


'am{kh/2y 


^ J-00 


kh/2 



coa{kx)dk. 



(5) 



These give progressively better approximations to the Gaussian at higher n, both by increasing the radius 
of compact support and by increasing smoothness, since each function Af„ is continuous up to the {n — 2}th 
derivatives. Since we minimally require continuity in at least the first and second derivatives, the lowest 
order B-spline useful for SPH is the M4 (cubic) spline truncated at 2h: 



w{q) = a 



i(2-g)3-(l-g)3, 

i(2-9)^ 



0<g< 1; 
1 < g < 2; 
9>2, 



(6) 



where for convenience we use VF(|r — r'|, h) = -^w{q), where q = \r — r'\/h and cr is a normalisation constant 
given by CT = [2/3, 10/(77r), I/tt] in [1,2,3] dimensions. Next are the M5 quartic, truncated at 2.5ft.: 



w{q) 



1) 'HI 



10 (i 



0<g< 
I <9< i; 
l<'z<i; 

<1>1 



(7) 



with normalisation a — [1/24, 96/11997r, l/207r], and the Mg quintic, truncated at 3/i: 



w{q) 



(3 




- 6(2 


~qf + Vo{l~q)\ 


< g < 1; 


(3 


-qf 


- 6(2 




1 < g < 2; 


(3 








2 < g < 3; 


0. 








'Z>3, 



(8) 
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with normalisation cr = [1/120, 7/4787r, l/1207r] (e.g. Morris 1996b). These kernel functions and their first 



and second derivatives are shown for comparison with the Gaussian in Fig. [2j 

One important aspect to draw from our discussion so far is the clear meaning attached to the smoothing 
length h as specifying the fall-off of the kernel weighting with respect to the particle separation. In particular, 
it is clear that referring to the "number of neighbours" does not have any meaning per se for Gaussian 
and Gaussian-like kernels: For the Gaussian the number of neighbours is in principle infinite, but there 
nevertheless exists a well-defined smoothing scale, h. The higher order B-splincs (Fig. [2| also demonstrate 
that it is possible to change the "neighbour number" - by progressing to higher n in the series - without 
changing the smoothing length. It is a widely-propagated myth that one can achieve formal convergence 
in SPH by "increasing the number of neighbours" (e.g. more slowly than the total number of particles). 
However, there are very important differences between simply "stretching" the cubic spline to accommodate 
a larger neighbour number - which amounts to changing the ratio of h to particle spacing - and using a 
kernel that has a larger radius of compact support but retains the same h. That is, in no sense is the SPH 
density estimate (our "approach 3") the same as approach 2 shown in Fig. [Ija. We will return to this point 
later. 

There are obviously other kernels, and other families of kernels, that satisfy the above properties (e.g. 
Dehnen 2001). However, more detailed investigations into kernels (e.g. Fulk and Quinn, 1996) tend to 



merely confirm the points made above - namely that Bell-shaped, symmetric, monotonic kernels provide 
the best density estimates. We will examine the formal errors in the kernel density estimate shortly, but 
first we turn to the issue of setting the smoothing length, h. 



2.4- Setting the smoothing length 



Early SPH simulations (e.g. Gingold and Monaghan 1977) simply employed a spatially constant resolu- 



tion length h, though one which was allowed to change as a function of time according to the densest part of 
a calculatiorj^ However, as is evident from Fig. [l] it is clearly desirable to resolve both clustered and sparse 
regions evenly - that is, with a roughly constant ratio of h to the mean local particle separation. Thus, a 
natural choice for setting the smoothing length is to relate to the local number density of particles, i.e.. 



h{r) (^n{r)-^/'^; 



ir) = Y,W[r-n,hir)] 

b 



(9) 



For equal mass particles, this is equivalent to making h proportional to the density itself (since 1/n = m/p). 
Since in turn density is itself a function of smoothing length, this leads to the idea of an iterative summation 
to simultaneously obtain the (mutually dependent) p(r) and h{r) ( Springel and Hernquist[ 2002 Monaghan 
2002 Price and Monaghanl 20071. Computed at the location of particle a, we have a set of two simultaneous 



equations 



P(l"a) 



b 



mbW{ra - n, ha)] 



h{ra) = 77 



Pa 



1/d 



(10) 



where 77 is a parameter specifying the smoothing length in units of the mean particle spacing (m/p)^/'*. These 
two equations can be solved simultaneously using standard root-finding methods such as Newton-Raphson 
or Bisection and most "modern" SPH codes employ such a procedure (for reasons that will become clear). 
Note that enforcing the relation given in ( 10 ) is approximately equivalent to keeping the "mass inside the 
smoothing sphere" constant (Springel and Hernquist 2002), since for example in three dimensions 



Ml, = / pdV. 



(11) 



^Similar to the spatially fixed but time-evolved gravitational softening lengths still employed in many cosmological simula- 
tions. 
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where Rkem is the kernel radius {2h for the cubic spline), so Mtot = const implies h'^ p = const. Since for 
equal mass particles Mtot = "m-Nneigh, this also means that the number of neighbours should be approx- 
imately constant if the relationship in ( 10 ) (or for unequal mass particles, Eq. [9]) is enforced. Indeed a 
"number of neighbours" parameter can be used in place of the parameter rj, using 



neigh,lD 



Km 



Nnetgh.2D = T^iCvY'i 



Nneigh,3D = T^T^iCvY 



(12) 



where ^ is the compact support radius in units of h (i.e., ( — 2 for the cubic spline). However, this is 
problematic for several reasons. Firstly it gives the dangerous impression that Nneigh is a free parameter 
unrelated to h, whereas changing Nneigh explicitly changes h ~ more specifically, the ratio of h to particle 
spacing (i.e., 77) - and corresponds to "stretching" the cubic spline as discussed above. Secondly, whereas 
r] carries the same meaning in 1, 2 and 3 dimensions, the Nneigh parameter changes, making it difficult 



neigh 

to relate the results of one and two dimensional test problems to three dimensional simulations 

- p (or h 



Nneigh is oftcn uscd as an integer parameter, whereas it is clear from ( 10 1 that the h 



Thirdly, 
n) iterations 

can be performed to arbitrary accuracy (that is, to fractional neighbour numbers) - which is also necessary 
if one is to assume that the relationship is differentiable. Finally, Nneigh is only related to the true number 
of neighbours so long as the (number) density of particles within the smoothing sphere is approximately 
constant (that is, so far as the integral in Eq. 11 can be approximated by |7ri?'^p). So at best a Nneigh 
parameter only characterises the mean neighbour number - and there can be strong fluctuations about this 
mean, for example in strong density gradients. 

Earlier adaptive SPH i mplement at io ns employed density estimates involving either an average smoothing 
length h — \ [hg + ht) (e.g. Benz 1990 1 or an average of the smoothing kernels Wab = \ \yVab(ha) + Wab{hb)] 



(e.g. Hernquist and Katz{ 1989p 



iowever, this inevitably leads to heuristic methods for setting the smooth- 
ing length itself - for example by evolving the time derivative of Eq. [TO] with "corrections" to try to keep 
the neighbour number approximately constant (e.g. Benz 1990 1 or simply enforcing a constant neighbour 
number either approximately or exactly (e.g. Hernquist and Katz 1989). It also considerably complicated 



attempts to incorporate derivatives of the smoothing length 
servation - into the equations of motion (Nelson and Papaloizou 1994). 



necessary for exact energy and entropy con- 
By contrast, the mathematical 



meaning of Eqs. ( 10 ) is clear and it is straightforward to take derivatives involving the smoothing length. 



Finally, the density estimate computed via ( 10 ) is time-independent, depending only on particle positions 



and masses and thus explicitly answering our original question of how to compute a density field from point 
mass particles. This also means it has wide applicability to many other problems beyond SPH - for example 



Price and Federrath (20101 use it to construct a density field from Lagrangian tracer particles in grid-based 



simulations of supersonic turbulence; and it forms the basis of the adaptive gravitational force softening 



method introduced by Price and Monaghan (2007). 



2.5. Errors in the density estimate 

The formal errors in the density estimate may be determined by writing the density summation as an 
integral - that is, assuming m = pdV and that the summation is well sampled, giving 



(p(r)) = / p{v')W(Y^v',h)dV'. 



(13) 



where (..) refers to a smoothed estimate. Expanding p(r') in a Taylor series about r, we have 



(p(r))=p(r) j W{r-r',h)dV'+\7p{r)-J{r'-r)W{r-r\h)dV'+\7°'\7^p{r) J Sr°' Sr^ W {r~r' ,h)dV'+0{h^) 

(14) 

so that if the normalisation condition ^ is satisfied and a symmetric kernel W{r — r', h) = W{r' — r, h) 
is used such that the odd error terms vanish, the error in the density interpolant is 0{h?). In principle it 
is also possible to construct kernels such that the second moment is also zero, resulting in errors of 0{h'^) 



(see Monaghan 1985). The disadvantage of such kernels is that the kernel function becomes negative in 
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some part of the domain, resulting in a potentially negative density evaluation. Achieving such higher 
order in practice also requires that the kernel is extremely well sampled, leading to substantial additional 
cost requirements. One possibility given the iterations necessary to solve ( 10 1 would be to automatically 
switch between high order and low order kernels during the iterations (e.g. if a negative density occurs), 
thus leading to high order interpolation in smooth regions but a low order interpolation where the density 



changes rapidly. The errors in the discrete version are discussed further in Sec. 4.3 



2.2 are not the only 



2.6. Alternatives to the SPH density estimate 

Finally, it should be noted that the three general approaches described in Sec. 
methods that can be employed for estimating the density. A fourth alternative that has received recent 
attention involves the use of Delaunay or Voronoi tessellation - the former proposed by Pelupessy et al. 
(2003) and the latter developed into a full hydrodynamics scheme by Serrano et al. (2005) and Hefi and 



Springel (20101. These are promising approaches that in principle can offer all the same advantages as SPH 
in terms of exact conservation - since it can be derived similarly from a Hamiltonian formulation - but with 
an improved density estimate and an exact partition of unity. 



3. From density to equations of motion 

The reader at this point may wonder why we have spent so long discussing nothing else except the density 
estimate. The reason is that this is the only real freedom one has if one wishes to obtain a fully conservative 
SPH algorithm, at least in the absence of dissipative terms. This is because the rest of the SPH algorithm 
can be derived entirely from the density estimate. 



3.1. The discrete Lagrangian 

The derivation starts with the discrete Lagrangian. As is usual, the Lagrangian is simply given by 

L = T~V, (15) 

where T and V are the kinetic and potential (in this case, thermal) energies respectively. For a system of 
point masses with velocity -v = dr/dt and internal energy per unit mass u, we have 



E 



TOfc 



-vl - Ub{pb, Sb) 



(16) 



where in general the thermal energy u can be specified as a function of the thermodynamic variables p 



and s (the density and entropy, respectively). Although (16) can be considered as a discrete version of the 



continuum Lagrangian for hydrodynamics (e.g. Eckart 1960 Salmon 1988 Morrison 1998 ) 



L 



[pv'^ - pu{p, s)] dV, 



(17) 



one is free to consider the discrete Hamiltonian system, it's associated symmetries and equations of motion 
directly - that is, without explicit reference to the continuum system. In other words the Hamiltonian 
properties are directly present in the discrete system and the motions will be constrained to obey the 
symmetries and conservation properties of the discrete Lagrangian. 

3.2. Least action principle and the Euler- Lagrange equations 

The equations of motion for such a system can be derived from the principle of least action, that is 
minimising the action 

S= f Ldt, (18) 
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such that SS ~ J 6Ldt ~ 0, where ^ is a variation with respect to a smaU change in the particle coordinates 
dr. Assuming that the Lagrangian can be written as a differentiable function of the particle positions r and 
velocities v, we have 



Integrating by parts, using the fact that 6v — d{5r)/dt, where d/dt = d/dt + v ■ W gives 

'dL 



SS = 



d 
dt 



dL 
dv 



dL 
dr 



Srydt 



dv 



■Sr 



0. 



(19) 



(20) 



to 



So if we assume that the variation vanishes at the start and end times, then since the variation Sr is arbitrary, 
the equations of motion are given by the Euler-Lagrange equations, here taken with respect to particle a: 



d_ / dL 
dt \dva 



dL 
<9r. 



= 0. 



(21) 



We have somewhat laboured the point here because it is important to understand the assumptions we have 
made by employing the Euler-Lagrange equations to derive the equations of motion. The first is that in 
using Eq. [21] we are not explicitly considering the discreteness of the time integral. So when we refer below 
to "exact conservation" (e.g. of energy and momentum) we mean "solely governed by errors in the time 
integration scheme"]^ 

The other, more critical, assumption we have made in employing (21 1 is that the Lagrangian is dif- 
ferentiable. This means that we have explicitly excluded the possibility of discontinuous solutions to the 
equations of motion. What this means in practice is that any discontinuities present in the system (for 
example in the initial conditions) require careful treatment - for example by adding dissipative terms that 
smooth discontinuities to a resolvable scale (i.e., a few h) such that they can be treated as no longer dis- 
continuous. We give some practical examples of this in Sec. |6.3| A be tter way would be to account for the 
neglected surface integral terms directly in the Lagrangian (e.g. Kats 2001 1, though it is not clear how one 
would go about doing so in an SPH context. 

3.3. Equations of motion 



All that remains in order to derive the equations of motion is to compute the derivatives in (21 



by 

writing the terms in the Lagrangian as a function of the particle coordinates and velocities. From ( 16 1 we 
have 

dpb 



dL 

dVa 



dL 
dra 



dub 
dpb 



dra' 



(22) 



the latter since m is a function of p and s, and we assume that the entropy s is constant (i.e., no dissipation). 
Note that the former gives the canonical momenta of the system (p = dL/dv). This step is straightforward 
for hydrodynamics, but it can also be used to derive the conservative momentum variable in the case of more 
complicated physics - for example in relativistic SPH (e.g. Monaghan and Price[ 2001 Rosswog 20091. 



3.3. L Thermodynamics of the fluid 

From the first law of thermodynamics we have 



dU = TdS - PdV, 



(23) 



where 5Q = TdS is the heat added to the system (per unit volume) and 5W = PdV is the work done by 
expansion and compression of the fluid. We do not compute the volume directly in SPH, but instead we 



''it should be noted that it is quite possible to also derive the time integration scheme from a Lagrangian — for example, 
[Monaghan] ( |2005| gives the appropriate Lagrangian for the symplectic and time-reversible leapfrog scheme. 
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can use the volume estimate given hy V = m,/ p and thus the change in volume given by dV ~ —m/p^dp. 
Using quantities per unit mass instead of per unit volume (i.e., du instead of dU), we have 



P 

du = Tds H — ^dp, 



(24) 



such that at constant entropy, the change in thermal energy is given by 



dub 
dpb 



P 



P 



(25) 



3.3.2. The density gradient 

So far we have not made any explicit reference to SPH or kernel interpolation. This arises because of the 



spatial derivative of the density (Eg. |22[), that we obtain by differentiating our density estimate (Eq. 10). 
Noting that we are taking the gradient of the density estimate at particle h with respect to the coordinates 
of particle a, this is given by 

dpb 1 ^ dWbc{hb) 
— — > rr 



dVa fib 



dVa 



{Sba — ^ca) , 



(26) 



where Wtdhb) = W{ri, — r^, hh), S^a is a Dirac delta function referring to the particle indices and we have 
assumed that the smoothing length is itself a function of density [i.e., h — h{p)], giving a term accounting 
for the gradient of the smoothing length given by 



dha v-^ dWab{ha) 



dpa 



dha 



where for the standard h — p relationship (Eq. 10 1 the derivative is given by 



dh 
dp 



h 

pd' 



(27) 



(28) 



where d is the number of spatial dimensions. 
3.3.3. Equations of motion 



Using (261 and (25) in (22) we have 



dL 

dVa 



Pb 



E 



dWUhb) . . 

rUc T. [Oba ~ Oca) : 



(29) 



which, upon simplification, gives the equations of motion from the Euler-Lagrange equations in the form 

Pa dWabiha) , Pb dW abijlb) 



dt 



b 



VtaPl dVa 



^bPl dVa 



(30) 



For a constant smoothing length these equations simplify to the standard SPH expression (Monaghan 1992 1 

'Pa , Pb 



dt 



J2mb 

b 



pI pI 



^aWab 



(31) 
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3.3.4- Conservation properties 

Although we cannot yet provide interpretation to the equations of motion so derived, we can note a 
number of interesting properties. The first is that the total linear momentum is conserved exactly, since 

J2 '"'^^'^ = J2''^^^ = -J2J2 ( 5 + 5 ) ^^^'^b = 0, (32) 

a a a b \Pa PbJ 

where the double summation is zero because of the antisymmetry in the kernel gradient]^ (the reader may 



verify that this is also true for Eq. 30). Secondly, the total angular momentum is also exactly conserved, 
since 

^TaX maVa = ^ma (^Ta X = ^ X! X! "^aTOft + ^ j Ta X (r,, - n)Fab = 0, (33) 



where for convenience we have written the gradient of the kernel in the form VaWab = ^abFab, and the last 
term is zero again because of the antisymmetry in the double summation [since (r^ x r;,) — — (rf, x r^)]. 
The above conservation properties follow directly from the symmetries in the original Lagrangian and 



(by extension) the SPH density estimate (Eq. 10 1 - linear momentum conservation because the Lagrangian 
and density estimate are invariant to translations, and angular momentum conservation because they are 
invariant to rotations of the particle coordinates. This is important in thinking about possible modifications 
to the SPH scheme (for example using non-spherical kernels would result in the non-conservation of angular 
momentum because the density estimate would no longer be invariant to rotations). 

Finally, we note that although the equations of motion depend only on the relative positions of the 



particles, they do depend on the absolute value of the pressure. That is, Eqs. ([30| and (31 ) contain a force 
between the particles that is non-zero even when the pressure is constant. We discuss the importance of - 
and problems associated with - this 'spurious' force in Sec. |5] 

3.4. Energy equation 

The remaining part of the (dissipationless) SPH algorithm - the energy equation - can also be derived 
from the Hamiltonian dynamics. Here we have the choice of evolving either the thermal energy u, the 
total specific energy e = -I- u or an entropy variable K = P/ . Equations for each of these can be 
derived, as given below. It is important to note, however, that - provided the equations are derived from 
the Lagrangian formulation - there is no difference in SPH between evolving any of these variables, except 
due to the timestepping algorithm. This is rather different to the situation in an Eulerian code where the 
finite differencing of the advection terms mean that writing the equations in conservative form (i.e., using 
e) differs more substantially from evolving u or K . 

3.4.1. Internal energy 



The evolution equation for the internal energy (in the absence of dissipation) , from Eq. 25 is given by 

dUa Pa dpa 



dt dt 



(34) 



Taking the time derivative of the density sum (Eq. 10), we obtain an evolution equation for u of the form 



^-7^5]mb(v,-Vb).V„W^„fc(/ia). (35) 



°This can be easily seen by swapping the summation indices a and b in the double sum and adding half of the original term 
to half of the rearranged term, giving zero. 
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3.4-2. Total energy 

The conserved (total) energy is found from the Lagrangian via the Hamiltonian 



H 



dL 



L 



1 



(36) 



which is simply the total energy of the SPH particles, i?, since the Lagrangian does not explicitly depend 
on the time. Taking the (Lagrangian) time derivative of (36), we have 

dE 



dt 



dUg 

dt 



(37) 



Substituting ( 30 ) and ( 35 1 and rearranging we find 



dE 



\ ^ dCa \ ^ \ ^ 



a 



a 



= 0. 



(38) 



This equation shows that the total energy is also exactly conserved by the SPH scheme (where the double 
sum is zero again because of the antisymmetry with respect to the particle index, similar to the conservation 
of linear momentum discussed above) . The conservation of total energy is a consequence of the symmetry 
of the Lagrangian ( 16 ) with respect to time as well as invariance under time translations. Eq. 38 also shows 
that the dissipationless evolution equation for the specific energy e is given by 

Pn _ Pb 



deg 
dt 



^mb 

b 



Wabiha) + T^V, • VaWabih) 



(39) 



3.4.3. Entropy 

For the specific case of an ideal gas equation of state, where 

P^K{s)p\ 



(40) 



it is possible to use the function K{s) as the evolved variable (Springel and Hernquist 20021, where the 
evolution of K is given by 



dK 

lit 



1 ■ 



1 



,7-1 



P 



du 
'dt 



The thermal energy is then evaluated using 



P dp 
'fPdt 



K 

7^ 



7 ■ 



1 



du\ 
dt) 



dis 



,7-1 



(41) 



(42) 



Since dK/dt = in the absence of dissipation, using K has the advantage that the evolution is independent of 
the time-integration algorithm. The disadvantage is that it is more difficult to apply to non-ideal equations 
of state. This is sometimes referred to as the 'entropy-conserving' form of SPH (after Springel and Hernquist[ 
2002 ) - which is somewhat misleading since the entropy per particle is also exactly conserved if (|35[) or (|39| 



are used provided the smoothing length gradient terms are correctly accounted for (i.e., du/dt — Pjfrdp/dt ~ 
0), apart from minor differences arising from the timestepping scheme. So the term 'entropy-conserving' more 
correctly refers to the correct accounting of smoothing length gradient terms and a consistent formulation 
of the energy equation than whether or not an entropy variable is evolved. 
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3. 5. Summary 

In summary, our full system of equations for p, v and u is given by 



Pa 

dt 

dug 
~dt 



^■mbW{va~rb,ha)] h = h{p), 

b 



^aPl 



^bp't 



P, 



^aPl 



^ ^ mb(vQ - Vb) • ^aWab{ha), 



(43) 
(44) 
(45) 



where in place of (45) we could equivalently use cither (39 1 or (41). The reader will note that so far we have 
not even mentioned the continuum equations of hydrodynamics - we have merely specified the physics that 



goes into the Lagrangian (Eq. 16 1, the thermodynamics of the fluid (Eq. 24) and the manner in which the 
density is calculated (Eq. 10) and with these have directly derived the discrete equations of the Hamiltonian 



system. In order to interpret them we need to know how to translate our SPH equations (43|-(45) into their 
continuum equivalents. 



Before doing so, it is worth recapping briefly the assumptions we have made in arriving at (43)-(45l 



from the Lagrangian (Eq. 16). These are 



i) That the time integration and thus time derivatives, d/dt, are computed exactly (though this assumption 
can in principle be relaxed); 

ii) That the Lagrangian, and by implication the density and thermal energies are differentiable; 

iii) That there is no change in entropy, such that the first law of thermodynamics du = PdV is satisfied, 
and that the change in particle volume is given by dV = —m/p^dp. 

The second and third assumptions in particular come into play in dealing with shocks and other kinds of 



discontinuities, which will we discuss further in Sec 6.3 



3.6. Alternative formulations 

Within the constraints of a Hamiltonian SPH formulation, it is clear that there is only a very limited 
freedom to change the algorithm without breaking some of the conservation properties of the scheme. 
So there are only two basic ways to change the (dissipationless) algorithm that are consistent with the 
Hamiltonian approach (other than changing the first law of thermodynamics): i) change the way the density 
is calculated or ii) introduce additional physical terms, and associated constraints, into the Lagrangian. 
Examples of the former are the consistent formulation of variable smoothing length terms - requiring the 



latter include consistent derivations of relativistic SPH (Monaghan and Price 2001 Rosswog 2009 ), adaptive 



gravitational force softening (Price and Monaghan 



2002[ |2004| ) and MHD ( ,Price and Monaghan, |2004b 



iterative solution of h and p in the density sum (Eq. 10 1 - by Springel and Hernquist (2002) and Monaghan 
(2002), and an incorporation of boundary correction forces ( Kulasegaram et al.[ 2004). Examples of the 



2007), sub-resolution turbulence models (Monaghan 



Price 2010) 



4. Kernel interpolation theory and SPH derivatives 

The usual way of introducing SPH is to start with a formal discussion of kernel interpolation theory. 
We have taken a rather different approach in this review and will introduce this theory primarily only to 
interpret the equations that we have derived from the discrete Lagrangian, and also as a way of discussing 
how to go about introducing additional physics. However, we will not use the linear error properties of the 
interpolation scheme to define the method - apart from our construction of the density estimate discussed 
above. The reason is that focussing on linear errors over the Hamiltonian properties of SPH misses some 
of the subtle but important non-linear behaviour that makes SPH work in practice, which we will come to 
discuss. However, first let us proceed: 
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.1. Kernel interpolation: the basics 
Kernel interpolation theory starts with the identity 



^(r) = j A{Y')6{r - r')dr', 



(46) 



where A is an arbitrary scalar variable and 5 refers to the Dirac delta function. This integral is then 
approximated by replacing the delta function with a smoothing kernel W with finite width i.e., 



where W has the property 



A(r) = J A{r')W{r~r',h)dr' + 0{h^). 



\imW{r -r',h) =(5(r-r'), 



(47) 



(48) 



and normalisation Jy WdV' = 1, as we have already discussed in Sec. |2.2l Finally the integral interpolant 

tides) b 



(Eq. 47 1 is discretised onto a finite set of interpolation points (the particles) by replacing the integral by a 
summation and the mass element pdV with the particle mass m, i.e. 

{A{v)) = j ^W{v~v\h)p{v')dv\ 



A, 



mb — W{r~~rb,h), 



(49) 



This 'summation interpolant' is the basis of all SPH formalisms. The reader will note that choosing A ~ p 
results in the SPH density estimate (Eq. 10). In this paper we have argued that the density estimate (10) is 
in some sense more fundamental than the summation interpolant ( 49 1 , since the equations of motion can be 
derived without reference to (49). On the other hand, the summation interpolant gives a general way of a 
interpolating a quantity at any point in space A{r) from quantities defined solely on the particles themselves 
(i.e., Ab, Pb, "rribl^ and in turn to a general way of formulating SPH equations. In particular, gradient terms 
may be straightforwardly calculated by taking the derivative of (49 1, giving 



VA(r) = 



dr 



^W{r~r',h)p{r')dr' + 0{h'), 
p(r') 

Ab 



ymb — VW{r-Tb,h). 



For vector quantities the expressions are similar, simply replacing A with A in (49), giving 

Ah, 



A(r) 
V • A(r) 
V X A(r) 
V^A^r) 



mb — W{r - Tb, h), 



E 



rub- 



Pb 



vw{v~n,h), 



nib— X Viy(r - r;,, h), 

b 



E 



7nb — Vm{r-rb,h). 

Pb 



(50) 
(51) 

(52) 
(53) 
(54) 
(55) 



The problem is that using these expressions 'as is' in general leads to quite poor gradient estimates, and we 
can do better by considering the errors in the above approximations. However, the basic interpolants given 
above give us a general way of interpreting SPH expressions such as those derived in Sec. |3] 



^Eq. |49| naturally also forms the basis for visualisation of SPH simulations, where one wishes to reconstruct the field in all 
of the spatial volume given quantities defined on particles: This is the approach implemented in SPLASH ( |Price[[2007[ l. 
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4-2. Interpretation of the Hamiltonian SPH equations 

We are now able to provide interpretation to our Hamiltonian-SPH equations (43H(45) derived in Sec.|3] 
using the basic identities (51)-(54). We begin with the density summation (43|. Taking the time derivative, 
we have 

^ = ^^mb(v,-Vb)-V„VKafc(/ia), (56) 
b 



dt n 

which for a constant smoothing length simplifies to 

dpa _ 

b 



dt 



= ^ nib (Vq - Vb) • WaWabih). 



Using (|53|) we can translate each of the terms according to 



dt 



Va • V Pb^aWab " V (PfeVfo) • V aWab ~ Va • V/3 - V • (pv) W -pa(V • v)a. 



Pb 



Pb 



(57) 



(58) 



So remarkably ( 57 1 (and by inference, 56 ) - which we obtained simply by taking the time derivative of the 
density sum - represents a (particular) SPH discretisation of the continuity equation. Indeed, the density 
summation is therefore an exact, time-independent, solution to the SPH continuity equation. This should 
not be particularly surprising, since the continuity equation derives from the conservation of mass, which is 
self-evidently enforced on the particles since m is fixed. 



Our force term (44), assuming a constant h (Eq. 31 1 can be translated according to 

P. 



b 



Pr, Ph 



Pi 



Pi 



p- 



rVp-V 



VP 

P 



(59) 



where we have used the basic identity (51 ) to translate the two terms into their continuum equivalents. 



Finally, the thermal energy equation (45 1 can be translated, from (34) and (58), giving 

du P_ 
— = V • V. 

dt p 



(60) 



In other words, it is evident that the Lagrangian has indeed given us valid discretisations for the equations 
of hydrodynamics, and therefore that the equations of our Hamiltonian system (43|-(45) do indeed solve 
these equations in discrete form - a remarkable achievement given the relatively few assumptions that were 
made. Yet the discretisations we derived are clearly not the basic ones arising from kernel interpolation 
theory (51)-(54). In order to examine the errors in these discretisations we need to understand the errors 
arising from the basic gradient operators, and how more general derivative operators can be constructed. 

J,. 3. Errors 



The errors introduced by the approximation (47) are similar to those in the density estimate (Sec. 2.5). 
pand A{v') in a 1 

A(r) + (r'-r)"- 



That is, if we expand A{y') in a Taylor series about r (Benz 1990 Monaghan 1992), we find 

dA 



{A{v)) = 



l(r'-r)"(r'-r)^^+0((r-rr 



W{r -r',h)dr', (61) 



such that for symmetric W = W{\r — r'|, h) and normalised (/ WdV = 1) kernels we have 

1 d^A 



(Air)) = A{r) 



2 9r"9r'3 



J (r' - r)"(r' - rfwi\r - r'|, h)dr' + 0[{r' - r)^ 



(62) 



giving an interpolation that is second order accurate [0{h'^)] unless higher order kernels are used (see 
Sec. 2.5). However, the errors in the discrete version (Eq. [49| are not identical, since they depend on the 
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degree to which the discrete summations approximate the integrals - specifically on the degree to which the 



discrete normalisation conditions are satisfied. Following Price (2004) we can perform a similar analysis on 



the summation interpolant (49 here assumed to be computed on particle a) by expanding At in a Taylor 
series around r^, giving 

(Aa) = V mh^Wab = V —Wab + VA, • V ^(r^ - ra)Wab + Oih'). (63) 

b b p" b P" 

This shows that in practice, the interpolation will only truly be second order accurate if the conditions 

— Wab « 1; and > 

Pb 

b ^° b 



V^(r,-r,)iy,fa«0, 
^ Pb 



(64) 



hold. The degree to which this is true depends strongly on the particle distribution within the kernel radius, 
and the properties of the kernel when a finite number of neighbours are employed - in particular, the ratio 
of smoothing length to particle spacing Ap. Fig. |3] shows the first condition computed for the B-spline 
kernels and the Gaussian as a function of h/Ap in ID (solid/black lines), showing that in general the above 
conditions are maintained very well, provided the particles are regular. Thus, maintaining a regular particle 
arrangement, together with an appropriate choice of h/Ap, can be very important in obtaining accurate 
results with SPH. Used another way, the conditions (64) are essentially the criteria for a 'good density 
kernel' - that is, a good density kernel is one which satisfies these conditions well given the typical particle 
distributions encountered in SPH. Note that the second of these is much easier to satisfy than the first, since 
it requires only a reasonably symmetric particle distribution to be satisfied. In general reasonable fulfilment 



of the first condition amounts to the conditions l)-3) on the density kernel described in Sec. 2.2 



The errors resulting from the gradient interpolation ( 50 1 may be estimated in a similar manner by again 
expanding A{r') in a Taylor series about r, giving 



VA(r) 



^(r) + (r' - r)"^ + ^(r' - r)^(r' - r)'',,;^ + 0[(r - r')^] 



Q,,Q,, ' -..^ .,,jVW^(|r-r'|,/.)dr', 
WWdr' + ^ j{v'- vrVWdv' + 2 j (r' - r)^(r' - vyVWAv' + 0[(r' - r)^]. 



(65) 



where we have used the fact that / VWdr' = for even kernels, whilst the second term integrates to unity 
for even kernels satisfying the normalisation condition ([s]). The resulting errors in the integral interpolant 
for the gradient are therefore also of 0{h^). As previously, the errors in the discrete version (51 1 can be 
found by expanding Ab in a Taylor series around Tq, giving 



y^rrib — VaWab, 
b P' 

. m;, dAa 



Pb 



rub 
Pb 



{n-rarVaWab+0{h^) 



(66) 



where the summations represent SPH approximations to the integrals in the second line of ( 65 ) . So the 



gradient errors resulting from (51) would in principle be similarly governed by the extent to which these 

(67) 



discrete summations approximate the integrals, i.e., how well the conditions 



V —VaWab 



0; and V !I^(rfc - r,)"Vf 1^,^ « 5"^; 



Pb 



hold. The latter term is shown for the B-spline kernels by the long-dashed/red lines in Fig.js] The difference 
is that for gradients we can explicitly use the error terms to construct more accurate gradient operators, 
and we do this below. 
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Figure 3: Accuracy with which the normalisation conditions on the kernel, the kernel gradient and the kernel second derivative 
are computed with fixed h on a one-dimensional line of particles with the M4-M6 kernels (from left) compared to the Gaussian 
(right), as a function of the ratio of smoothing length to particle spacing h/Ap. Bell-shaped kernels with compact support 
lead to excellent density estimates (solid/black lines), reasonable gradient estimates (long dashed/red lines) but poor second 
derivative estimates (short dashed/green lines) when a finite number of neighbours are employed. 



4.4- First derivatives 

From (66) we immediately see that a straightforward improvement to the gradient estimate (51 1 can be 
obtained by a simple subtraction of the first error term (i.e., the term in (66 1 that is present even in the 
case of a constant function), giving (e.g. Monaghan, 1992) 



b 



Pb 



-^aWab, 



which, interpreting each term according to (51 1, is an SPH estimate of 

(VA)-A{V1). 



(69) 



Since the first error term in (66) is removed, the interpolation is exact for constant functions and indeed 
this is obvious from the form of ( 68 ) . The interpolation can be made exact for linear functions via a matrix 

«/3 _ "^fc , 



inversion of the second error term in (66), i.e., solving 



Pb 



b P" 



(70) 



where = d/dr^ . This normalisation is somewhat cumbersome in practice, since x is a matrix quantity, 
requiring considerable extra storage (in three dimensions this means storing 3x3 = 9 extra quantities for 
each particle) and also since calculation of this term requires prior knowledge of the density. 
A similar interpolant for the gradient follows by using 



« - [{V(pA)) - A{Wp)] = - V mfc(Afc - Aa)VaWab 

p Pa. 



(71) 



which again is exact for a constant A. Expanding Af, in a Taylor series, we see that in this case the 
interpolation of a linear function can be made exact by solving 



x''^^-Y.^biAb~Ag)V^Wgb, 



X 



Y.^biVb-VaT^f'Wab. 



(72) 
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which has some advantages over ( 70 1 in that it can be computed without prior knowledge of the density. 



However, the gradient operator we derived in the equations of motion (31) does not correspond to any 
of the above possibilities. This operator is given by 



-^(Vp) + (V' 



Expanding Af, in a Taylor series about r^, we have 



Pa'^mb 

b 



Pi 



Ah , 

Pb 



PaA^ ^ ™ J 1 + 1 ) VaWab + -JT^ Pa E 5 " rJ"V,W^,, + 0{h^ 
J \Pa Pb/ b Pb 



from which we see that for a constant function the error in (73 1 is governed by the extent to which 

(4 + 4 ) ^aWab - 0. 
b \Pa Pb J 



(73) 



(74) 



(75) 



Although a simple subtraction of the first term in ( 74 ) from the original expression ( 73 1 eliminates this error 



this would not give the form we derived in Sec. 3.3 Indeed, retaining the exact conservation of momentum 
requires that such error terms are not eliminated, the consequences of which we will discuss in Sec. [Sj 

^.5. Generalised first derivative operators 

Finally, an infinite variety of gradient operators - of two basic types - can be constructed by noting that 



VA = i [V(M) - ^V0] « V — $^ (A - Aa) VaWab, 
<t> "h" Pb (pa 



(76) 



and 



E 



rub 
Pb 



-Aa 



'^aWab. 



(77) 



where <f> is any arbitrary, differentiable scalar quantity defined on the particles. Indeed, for a given <j), the 
pair of operators defined by (76) and (77) can be shown to form a conjugate paiij^- and choosing one (e.g. 
for the density /thermal energy evolution) tends to lead to the other (e.g. in the equations of motion). For 
example, (68) and (71) correspond to using (j) — \ and <f) — p respectively in ( |76[ ), whilst ( |73[ ) correponds to 
using (p ~ p in (77) and arises in the equations of motion because it is the conjugate operator to (71 ) that 



arises in the density gradient (57) 



Various 'alternative' formulations of the SPH equations have been proposed that correspond to a par- 



ticular choice of (j) (e.g. 


Hernquist and Katz 


1989 Ritchie and Thomas 


2001 


Marri and White 


2003 


example 


Hernquist and Katz 


(1989 


) suggested using an acceleration equation of the form 



'dt 



J2mb 
b 



2^1™ ) ^aWab, 
PaPb 



(78) 



corresponding to the symmetric operator (Eq. 77) with = \[P j p. It can be readily shown that ensuring 
exact conservation of energy with such a formulation would require using the conjugate operator (i.e., 
Eq. 76 with = 



P I p) in the thermal energy equation (although it should be noted that simultaneous 



exact conservation of energy and entropy is not possible with any alternative formulation) . 



^The conjugate nature of the symmetric and antisymmetric SPH gradient operators was first noted by [Cummins and| 
iRudmanl ||1999 \ in the context of projection schemes for SPH. 
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4-6. First derivatives of vector quantities 



All of the above discussion applies also to vector derivatives, simply using (53) or (54) in place of (51 1 



and likewise resulting in two basic operators for each type of derivative, given by 
(V.A),« V"^$^(A,-A„).V„Ty„fc; (V x A>„ 

^ Pb (Pa 

and 

TUb ( (j)b 



^"^$^(A,-A„)xV,W^,,, (79) 

^ Pb (Pa 



"T^ Pb \(Pa 



06 



^aWab 



^ Pb \(Pa 



Ab X WaWab 



(80) 

where as previously is an arbitrary (difFerentiable) scalar quantity. For general vector derivatives written 
in tensor notation the corresponding expressions are given by 

mb 



^ Pb 0a ^ Pb \(Pa (t>b 



rUb <j>b 
Pb (tia 

These operators form the basic building blocks for formulating quite general SPH equations. Higher order 



operators can also be constructed for vector derivatives using matrix inversions, similar to (70) and (72). 



4. 7. Particle methods "the wrong way": an SPH formulation based on linear errors and exact derivatives 

Based on the discussion given in Sec. |4.3H4.6| a straightforward approach would be to simply go ahead and 
discretise the continuum equations of hydrodynamics (or magnetohydrodynamics) using the most accurate 



gradient estimates possible. For example employing (68) the hydrodynamic equations of motion dv/dt 
—VP/ p could be written in the form 



E^"h 



PaPb 



^aWab, 



(82) 



or any of the alternatives offered by choosing cj) appropriately in ( 76 ) . Indeed such a formulation was originally 
examined (and discarded) by Morris ( 1996a|b I but has been recently (re-)proposed by Abel ( 2010 ) (the latter 
author employing (j) — 1/p). We could even proceed to more accurate derivatives using Eqs. (70) or (72 1 
that are exact to linear order. Yet, these operators are clearly different from the equation of motion we 



derived from the Lagrangian (30), though on the basis of the linear error properties alone (Sec. 4.4) would 
seem to be a much better choice. On the other hand, it is clear that (82) does not exactly conserve linear 
(or angular) momentum, nor total energy. 

The key difference in approach is that the above analysis tells us about the linear errors, whereas the 
Hamiltonian formulation tells us about the non-linear properties of the system (i.e., symmetries, conservation 
and constraints on the behaviour of the global system). These turn out to be crucial for long term stability 
and accuracy, and we will look at what this means in practice in Sec. [S] That said, "linear errors do not lie" , 
so it will always be true that - provided the particle distribution is regular - formulations such as ( [82| ) or 
similar will be more accurate for linear or weakly non-linear problems (i.e., those run for a short time and/or 
not involving strong shocks) and with sufficient resolution can always be made to give accurate results. 

Ideally of course, one would want both exact derivatives and exact conservation. In SPH at least, it 
seems one cannot have both - and to my knowledge this is yet to be convincingly demonstrated by any 
particle method, though it is perhaps possible using tesselation schemes. 



5. Why a bad derivative leads to good derivatives: The importance of local conservation 

The paradox we face is that, whilst the Lagrangian derivation gave us valid discretisations of the equations 
of hydrodynamics, these are not the most accurate discretisations that are possible based on a linear error 
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analysis. Indeed, on the basis of the linear error properties, the gradient estimate derived for the acceleration 
equation would be a very poor choice of gradient operator, yet it is the only operator which respects all of 
the symmetry and conservation properties in the Lagrangian. 

To understand what these errors mean in practice we can perform a simple thought experiment (which 
we will also run as a numerical experiment): Consider a distribution of particles in a closed (e.g. periodic) 
box at constant pressure. In particular, we will consider a uniform random particle distribution. This is 
the worst-case error scenario, since errors in the interpolation will essentially be Monte-Carlo (cx 1/\/N). 
Now consider what will happen if use one of the more accurate gradient estimators, for example (82). 
Clearly, since the pressure is constant, there will be no force and hence no particle motion. That is, the 
force vanishes when the pressure is constant regardless of the particle distribution. If instead we use the 
Hamiltonian formulation, then the pairwise force between particles is given by 



-rriamb 



P P 



Pi 



Pi 



Fab (fa 

—ve 



(83) 



where we have written the kernel gradient using VaWah = ^abFab- Since for positive definite kernels (Fig.[2| 
the kernel derivative function is always negative, the overall result ~ assuming positive pressure - is a positive 
(that is, repulsive) force between the particles directed along their line of sight. This force will only vanish 
when the error term (75) becomes zero - that is, when the particles are regular. What this means is that the 



particles care about their (bad) arrangement and will rearrange themselves until the condition (75) holds, 
corresponding to a particle distribution that is locally isotropic and regular. This state corresponds to the 



minimum energy state of the Hamiltonian system - i.e., that which minimises the action (181, and in this 
state the symmetric gradient estimate ( |73[ ) is computed with good accuracy (since by definition the motions 
act to minimise the errors in this estimate). 



Monaghan (2005) gives some specific examples of this. 



In other words, the Hamiltonian SPH formulation - or more generally any formulation which respects 
local conservation of momentum between particle pairs - contains an intrinsic "re-meshing" procedure and 
the particles are constrained to remain locally ordered at all times. With an 'accurate' but non-conservative 
pressure estimate the particles are insensitive to their arrangement and can become arbitrarily randomised 
in the course of a simulation (according to the streamlines of the flow). Thus, methods based on an 'exact 
derivatives' approach (e.g. Dilts 1999 Maron and Howes 20031 inevitably require the addition of explicit 
(and ad-hoc) re-meshing procedures, as the interpolation errors on a randomised particle distribution will be 
terrible regardless of the order of the interpolation scheme. So in practise the 'accurate' gradient estimate 
will give much less accurate results. This is the paradox of SPH: One deliberately chooses a bad gradient 
estimate (in the linear errors) in order to obtain a good gradient estimate (because the particles stay regular). 

One of the implications of this "intrinsic remeshing" is that not all initial conditions represent stable 
configurations for the particles. In particular, the cubic lattice - though simple - does not represent a 
minimum energy state and is only quasi-stable for certain ratios of the smoothing length to the particle 
separation (Morris 1996b B0rve et al. 2004). In general, given a small perturbation or sufficient time, the 



cubic lattice will transition to the more isotropic hexagonal-close-packed lattice arrangement and will do so 
almost immediately if /i/Ap is small (e.g. rj < 1.1). 

5.1. Example 1: Settling of a random particle distribution 

Our first numerical example (using ndspmhd) demonstrates this 'settling' in practice. The setup is a 
two dimensional domain x,y d [0, 1] with periodic boundary conditions. The particles are given an initially 
uniform thermal energy, density is calculated according to the sum and the pressure is determined using 
P = (7 — l)pu with 7 — 5/3. The thermal energy is set to give a sound speed = 1, giving u — 0.9. 
The end result should therefore be a uniform pressure equilibrium. The particle distributions after 0.5 (top 
row) and 10 (bottom row) sound crossing times are shown in Fig. |4]using the Hamiltonian SPH formulation 
(43|-(45) (left and centre columns, without and with artificial viscosity respectively) and the equations of 
motion computed using a 'relative pressure' formulation (specifically, Eq. 76 with = l/p as employed 
by Abel 2010| ) (right column). With a locally conservative formulation the random initial configuration 
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Figure 4: Settling of an initially random particle distribution due to pairwise conservation of momentum in the pressure 
gradient. The left and centre columns show the results with the standard (Hamiltonian) SPH method, both without (left) and 
with (centre) artificial viscosity, after 0.5 (top) and 10 (bottom) sound crossing times. The right column shows the results when 
a "relative pressure" formulation is adopted. With a momentum-conserving force the particles are sensitive to their arrangement 
and will regularise accordingly (left and centre columns), whereas "more accurate" but non-conserving formulations (right) 
compute gradients that are insensitive to the particle arrangement and thus require explicit re-meshing procedures. Note that 
although the application of artificial viscosity helps the settling to proceed faster (centre), it is primarily a pressure-driven 
effect and occurs even if no viscosity is applied (left). 



settles rapidly into a regular particle distribution (left and centre columns), leading in turn to good gradient 
estimates. This settling occurs even in the absence of an artificial viscosity term (left panels), though adding 
viscosity does help speed the settling process (centre panels). However, using a "more accurate" but non- 
conservative gradient estimate there is no regularisation of the particle distribution, leading to poor gradient 
estimates due to the random nature of the particle distribution. The total energy is also conserved exactly 
by the SPH formulations, whilst in the relative pressure formulation the total energy grows exponentially. 



5.2. Example 2: A 2D shock tube 

The other "classic" example of particle settling is the behaviour of SPH particles in a multidimensional 
shock tube problem, where there is a ID compression of the particle distribution (e.g. along the x axis). Since 
the shock induces a highly anisotropic compression - and thus a highly non-preferred particle arrangement 
- the mutual repulsion of SPH particles will eventually produce a post-shock "remeshing" of the particle 
distribution, involving transverse motions of the particles. An example is given in Fig.[5j showing the particle 
distribution at t = 0.1 in a two dimensional Sod shock tube problem (described further in Sec. 6.3.4) in 
which the particles were initially placed on two hexagonal close packed lattices upstream and downstream 
of the shock (initially placed at the origin). The particles can be seen to "break" (at x ~ 0.14) from 
the highly anisotropic compression-induced 'lines' at 0.14 < a: < 0.19, leading to a more isotropic particle 
distribution further downstream [x < 0.12). This example also illustrates the fact that one inevitably has 
some motions at the resolution scale that are not related to the physical problem, but related to the implicit 
"regularisation" of the particle distribution present in locally conservative SPH formulations. 
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Figure 5: Particle settling in a two dimensional shock tube problem. The particles are initially arranged on hexagonal close 
packed lattices either side of the shock (top panel). As the shock propagates (bottom panel, showing t = 0.1) it induces a 
one dimensional compression in the particles and thus a highly anisotropic particle arrangement, which "remeshes" to a more 
isotropic arrangement downstream from the shock, involving small motions of particles in the y— direction. 



5.3. Corollary: Negative pressures and the tensile instability 

The corollary of the above is that the particles require a positive pressure in order to remain ordered. If 
the net pressure (or stress) becomes negative, the net force between a particle pair will become attractive, 
causing a catastrophic numerical instability. For example, with a pressure gradient of the form 
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(84) 



the pairwise force will become negative when Pq > P^ and in this situation the particles clump together 



unphysically. This is known as the 'tensile instability' (Monaghan 2000) and occurs in SPH when a stress 



tensor is employed that can result in (physically) negative stresses. In particular, this is the case for MHD 
( Phillips and Monaghan 1985) and in elastic dynamics (Gray, Monaghan, and Swift 2001). The occurrence 



of the tensile instability was one of the main initial difficulties with the development of MHD in SPH and 
is discussed in detail in Sec. |8l 



5.4. The pairing instability: Why one cannot simply use 'more neighbours'. 

Another, more benign, instability in the particle distribution occurs with the cubic spline and other 
bell-shaped kernels depending on the ratio of smoothing length to particle spacing. This is due to the shape 
of the kernel gradient term for these kernels (see Fi g. [2| ), and is a consequence of the fact that these kernels 
are designed to give good density estimates (Sec. |2.2| , rather than necessarily being the best choice for 
calculating gradients. In particular, the kernel gradient in these kernels contains a maximum (negative) 
value at r/h ~ 2/3 and tends to zero at the origin (Fig. [2]). This characteristic is desirable for a good 
density estimate - as it means one is insensitive to a small change in the position of a near neighbour - but 
means that the mutual repulsive force tends to zero for neighbouring particles placed "within the hump" of 
the kernel gradient. The net result is that two particles spaced closer than the location of the "hump" in 
the gradient form a "pair" , eventually falling on top of each other. 

For the cubic and other B-spline kernels complete merging occurs when h > 1.5Ap (i.e., rj > 1.5 or 
> 100 Neighbours in 3D for the cubic spline), corresponding to the placement of the first neighbour "inside 
the hump". There is also an intermediate regime 1.225 ^ ?7 ^ 1-5 (62-100 Neighbours in 3D) where a 
close-packed or cubic lattice is unstable to pair formation, but where the pairs do not completely merge. 



These empirical regimes are confirmed by detailed stability analysis of the SPH equations in 2D (Morris 
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Figure 6: The pairing instability in action: The (2D) setup is similar to that shown in Fig. |4] except that the particles are 
initially placed on a close packed lattice and we use the M4 cubic spline kernel with a large ratio of smoothing length to particle 
spacing (here rj = 1.5, corresponding to ^ 28 and 100 neighbours in 2 and 3D respectively). After a few sound crossing times 
(t = 5, centre panel) particles form 'pairs' which proceed to merge into a locally hexagonal "glass-like" lattice arrangement with 
almost exactly half the resolution of the initial conditions (right panel, shown after 10 sound crossing times). Although fairly 
benign — and easily avoided by a sensible choice of jj - the pairing instability is the main reason one cannot simply "stretch" 
the cubic spline to large neighbour numbers to achieve convergence. Instead, one should use a kernel with a larger radius of 
compact support but the same ratio of smoothing length to particle spacing, such as the M5 or Mq splines. 



1996b[ B0rve et al. 2004) that explicitly show that instability occurs - though with small energies - for 



large h/Ap. 

Fig. [6] (and our example 3) shows the pairing instability in action: The setup is as for example 1 but with 
= 1.5 in the cubic spline kernel instead oi rj — 1.2 and with particles placed initially on a hexagonal close- 
packed lattice (an otherwise very stable configuration: left panel). After a few sound crossing times (centre 
panel) particles begin to form pairs, with these pairs eventually merging completely (right panel) to give a 
locally hexagonal "glass-like" configuration, but with exactly half the resolution of the initial conditions! 

Though fixes have been proposecj^ none are entirely satisfactory. However, the pairing instability, unlike 
the tensile instability, is quite benign. For example the density change associated with the transition in 
Fig. [6] is of order 1% for the cubic spline and 0.1% for the Mg quintic - but entails a factor-of-two loss in 
spatial resolution and is therefore a waste of computational resources. Furthermore, it can be easily avoided 
by a sensible choice of t] (we recommend rj = 1.2 for the B-spline kernels, corresponding to Nndgh — 57.9 
for the cubic spline in 3D). The pairing instability is the main reason one cannot simply "stretch" the cubic 
spline to large neighbour numbers in order to obtain convergence and demonstrates at least one good reason 
why rj (or Nndgh) should not be regarded as a free parameter in SPH simulations. 



^Thomas and Couchmaii| |1992[ | suggested modifying the gradient of the cubic spline kernel, using 

(-1 0<g<2/3; 
-3g+|q2, 2/3<g<l; 
f(2-<?>, 1<,<2; 
0. q>2. 

with W itself unchanged and a equal to the usual normalisation factor for the cubic spline (i.e., l/vr in 3D). That is, the "hump" 
is removed by simply making the kernel gradient constant within r/h< 2/3. Whilst it cures the pairing instability, one should 
be careful about employing such a gradient in practice since the kernel gradient | |85| l is no longer correctly normalised (i.e.. 
Eg. |67[ 3 no longer holds, even in the continuum limit) meaning that as the region within r/h < 2/3 is increasingly well sampled 
the numerical sound speed and other quantities will be systematically wrong. Though one could attempt to re-normalise the 
new gradient kernel, this results in a low weighting in the outer regions that in turn leads to poor gradient estimates. 

Similarly, whilst perhaps a satisfactory 'gradient kernel' could be derived without a pairing instability, in the derivation from 
a Lagrangian there is no freedom over the kernel gradient since it derives directly from the gradient of density — that is, if one 
separates the gradient kernel from the density kernel then either the total energy (fr om Eq. |37| or the entropy will no longer 
exactly be conserved (the latter if du/dt 7^ P/ p^dp/dt and thus dK/dt ^ in Eq.pll. 
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6. Second derivatives and dissipation terms in SPH and SPMHD 



6.1. The SPH Laplacian 

Our remaining "basic" issue regarding botlr SPH and SPMHD regards the formulation of second deriva- 



tive terms. As for first derivative terms we can start witli tire basic summation interpolant (49) and simply 
take derivatives analytically, e.g. 

{V^A)a^J2'^b^^l^-b. (86) 
b 



Pb 



Expanding Af, in a Taylor series about r^, we have 
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Pb 



Pb 



Pb 



(87) 

As previously, we can immediately improve the estimate by subtracting the first error term from both sides, 
giving an interpolant of the form 

Y—^^^- Aa)^l'^ab = V" A„ Y —S^'^'^lWab + J V" ^ ^ 1^ 5v'^ Sv^ ab + ■ ■ • , 
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which vanishes when A is constant. However, the accuracy of our second derivative estimate will depend on 



the remaining error terms in (881, corresponding to the normalisation conditions: 



such that 
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The problem is that the conditions ( 89 1 are very poorly satisfied using the second derivative of a compact 
bell-shaped kernel function (short-dashed/green lines in Fig. [2j showing the second term in 89), since the 
second derivative changes sign inside the kernel domain (for the cubic spline it is also discontinuous) and 
must therefore be extremely well sampled in each direction to give accurate results. Thus in practice we 
require a better functional form - a "second derivative kernel" . 

The actual formulation commonly employed derives from early work by Monaghan and first published 



by Brookshaw (1985), where the SPH Laplacian is written in the form 



V^Aa^2Y—{Aa~A, 



Pb 



Fab 
Vab\ 



(91) 



where we use the definition VaM^a& = i'abFab such that Fab is the scalar part of the kernel gradient term. 
Thus in effect we use the first derivative kernel function divided by the particle spacing to give a second 



derivative. Whilst the interpretation of (91 ) from the integral representation is in general rather complicated 



(see, e.g. Espafiol and Revenga 2003 Jubelgas et al. 2004 Price 2004 Monaghan 2005), it can be more 
easily understood as if we had simply used (90) with a new kernel Yab instead of Wab, defined according to 



y^Yab 



'^Fab 

Kb\ ' 



(92) 



The left figure in Fig. [t] shows Y" constructed in the above manner from the A/4 and Mq kernel gradients, 
and may be compared with the standard kernel second derivatives shown in Fig. [2] The right figure shows 
the normalisation condition (89 1 for these 'kernels'. Comparison with Fig. [3] shows that, indeed, the second 
derivative is much better estimated with kernel second derivative functions that are monotonically decreasing 



and positive, such as those constructed via (92) from the bell-shaped kernels 
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Figure 7: Second derivative kernel functions Y"(q) = —2w'(q)/q constructed from the first derivatives of the M4 cubic and Mg 
quintic kernel functions (left figure, compare to the stand ard second derivatives shown in Fig. [2]l, together with the accuracy 
with which the second derivative normalisation conditions (jsojl are satisfied for a fixed ratio oi h/Ap (compare with the standard 
second derivative functions shown in Fig. [sj . 



Understanding the 'BrookshawJ ( 1985 ) Laplacian as equivalent to (90 1 with an ahernative kernel also 
helps to interpret more complicated expressions. For example, it is quite straightforward to show that 



Y.—{Ka + Kb){Aa-Ab) 
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(93) 



by writing —2Fab/\rab\ = V^Fob and interpreting each term via (86 1. Cleary and Monaghan ( 1999 1 proposed 
an alternative average of the k terms given by 



ErUh 4:KaKb , . A \-^<^'> 

7 I \\^a — ^b)-, r 
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(94) 



which was formulated to give smooth derivatives when n is discontinuous. The above expression forms 



the basis of formulations of thermal conductivity in SPH (Cleary and Monaghan 1999), though has been 
similarly used to model a wide range of dissipative terms including viscosity (e.g. Cleary et al. 2002 Hu 



and Adams 2006), salt diffusion (Monaghan, Huppert, and Worster |2005) and in an astrophysical context 



for the treatment of radiation in the Flux-Limited Diffusion approximation (Whitehouse and Bate 2004 



Whitehouse et al. 2005). 



6.2. Vector second derivatives 

Second derivatives of vector quantities do not quite follow the same analogy as the Laplacian, because in 
general V*V^ W^ab involves a mix of the first and second derivatives of the dimensionless kernel function (see 



appendix A in Price 2010). Thus, the proof is more involved (see Espahol and Revenga 2003 Monaghan 



2005 for details), but the basic expressions for vector second derivatives are given by 
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where = d i.e., the number of spatial dimensions. A corollary of the above is that a second derivative 
computed purely along the line of sight between the particles (e.g. constructed so as to conserve angular 
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momentum) corresponds to 
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(97) 



An alternative approach is to calculate vector second derivatives by taking two first derivatives. Although 
more expensive, this has been the approach adopted in several formulations of physical viscosity, e.g. for 



SPH modeUing of accretion discs (Flebbe et al. 1994 Watkins et al. 1996) 



6.3. Artificial dissipation terms in SPH and SPMHD 
6.3.1. Interpretation of SPH artificial viscosity terms 

This brings us to the formulation and interpretation of artificial dissipation terms in SPH. The 'standard' 



formulation of artificial viscosity is given by (Monaghan 1992) 
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(98) 



where barred quantities correspond to an average, i.e., pab = {Pa + Pb)/'2', Cg is the sound speed, a and /? 
are dimensionless parameters (typically a = 1 and /3 = 2) and e ^ 0.01 is a small parameter to prevent 
divergences. This form was chosen because it is Galilean invariant, vanishes for rigid body rotation and 



conserves total linear and angular momentum (Monaghanj 1992 1. If we neglect the non-linear (/3) term, set 
e = and assume that Cs, h and p are approximately constant over the kernel radius, then using (97 1 we 
can directly translate this expression into the continuum form 
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(99) 



Comparison with the compressible Navier-Stokes equations shows that the artificial viscosity is therefore 
equivalent to a physical viscosity with shear and bulk coefhcients proportional to the resolution length, i.e.. 



(e.g. Artymowicz and Lubow 1994 Murray 1996 Monaghan 2005 Lodato and Price 20101 
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(100) 



Since for shock-capturing the viscosity is only applied when particles are approaching (vat, • Tab < 0), in a 
uniform shear flow - or accretion disc ~ the viscosity coefficients will be approximately half of these values. 

6.3.2. General formulation of dissipative terms in SPH and SPMHD 



A more general formulation of dissipative terms was proposed byjMonaghan ( 1997), based on an analogy 
with Riemann solvers and the need to formulate dissipative terms for ultra-relativistic shocks (Chow and 



Monaghan 1997 ) . The general principle is that dissipative terms in the conservative variables involves jumps 



in those variables (the "left" and "right" states of the Riemann problem), multiplied by eigenvalues that 
can be interpreted as signal velocities. 

For the equations of hydrodynamics the conservative variables are the density, specific momentum a nd 
energy (p, v and e — ^v"^ 



respectively), and the terms take the fornj^ (Monaghan 

stgi'Va - Vfe) • iab 
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(102) 



^Note that in SPH no dissipation is required in the density equation provided the density is comp ute d from a sum. This is 
because represents an integral form of the continuity equation and can be shown to differ from ( |56| by a surface integral 
term that is non-zero only at boundaries or equivalently, discontinuities in the flow. See [Price for more details. 
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Figure 8: Treatment of discontinuities in SPH: A ID Sod shock tube problem (left) and 2D Kelvin- Helmholtz instability (right), 
applying only artificial viscosity terms and with unsmoothed initial conditions. Whilst in the ID shock tube problem (left 
panel) the shock is smoothed over several particle spacings by the viscosity term, there are problems at the contact discontinuity 
causing a 'blip' in the pressure. This leads to a suppression of mixing across contact discontinuities in 2D (right panel) 





Figure 9: As in Fig. [8] but with artificial conductivity applied as well as artificial viscosity. This smoothes the contact 
discontinuity, removing the pressure 'blip' (left panel) and restoring the mixing in the 2D K-H problem (right panel). 



(vo • fafc)^ + ctuV^^gUa refers to an energy including only components along the line of 
sight joining the particles, Fab = [Fab{K) + Fab{hb)]/2 (equivalent to writing Fab = ^ab ■ ^aWab) and a and 
Uu are the dimensionless artificial viscosity and thermal conductivity parameters, respectively. The signal 
speed Vsig refers to the maximum (averaged) signal speed between a particle pair. For hydrodynamics we 
use 



2 I'^s.a 
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b ■ TafcJ 
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Vab • fab > 0. 



(103) 



Eq. (101 1, with the above Vsig, provides a standard artificial viscosity term similar to ([98|) (careful expansion 



shows they differ only by a factor of h/\rab\ and the no- longer-necessary e term). The dissipation term in 
the total energy also contains a term involving (ua — Ub) which acts to smooth jumps in the thermal energy 
(e.g. at a contact discontinuity). The interpretation of this term is clearer from the contribution to the 
thermal energy evolution, given by 
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Figure 10: The Sod shock tube in 2D. Using the cubic spline kernel (left panel), there is additional 'noise' in the 2D velocity 
field (compared to ID, Fig.[9]| due to the transverse "remeshing" motions of particles behind the shock front in multidimensions 
(see Fig.|5]|. However, this can be quite effectively minimised by using a smoother kernel (right panel, using the Ma quintic). 



where the second term, from ( 91 1, can be directly mterpreted as a V^m term. Note that the signal velocity v^^^ 



used in. the artificial conductivity does not have to be the same as that used for the viscosity. In particular, 

Pb\/Pab to equalise the pressure across contact discontinuities, 

|vaf, ■ fa&l (we adopt 



Price 
whilst 



( 2008 1 proposed using v^^g = \/\P^ 



Wadsley et al. (2008) proposed a conductivity term equivalent to using v'^^^ 
the former tor the tests shown in this paper). 



6.3.3. Switches for viscosity terms 

One of the key issues in practice is to ensure that sufficient dissipation is applied to discontinuities, but 
that such dissipation is effectively turned off in smooth parts of the flow by designing appropriate switches. 



Morris and Monaghan (1997) suggested allowing the parameter a to be individual to each particle, with an 
evolution equation of the form 



da 
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(105) 



where 5 is a source term that grows large at the discontinuity [e.g. S ~ max(0, — V • v) for shocks], r 
is the decay time, set such that a decays to over several smoothing lengths (typically cr = 0.1 and 

nin = 0.1). A similar switch can be employed for the thermal conductivity parameter a„, with [Price and 
Monaghan (2005) adopting a source term given by Su = O.lhV^u. More sophisticated switches for shock 
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6.3.4. Examples 4 cind 5: One and two dimensional shock tubes, and Kelvin- Helmholtz instabilities 

Two specific examples of the dissipative terms in practice are shown in Figs, [s] (applying only viscosity) 
and [9] (applying artificial viscosity and conductivity), showing the results of a one dimensional Sod shock 
tube problem (left subfigure) and a two dimensional Kelvin-Helmholtz (K-H) instability problem (right 
subfigurc). The ID shock tube is setup using a total of 450 particles in —0.5 < a; < 0.5 with conditions to 
the left of the origin given by [pL, Pl,vl] = [1,1,0] and to the right by [pr, Pr,vr] = [0.125,0.1,0] with 
7 = 5/3. Importantly, purely discontinuous initial conditions are employed so that the contact discontinuity 
is not already smoothed. The K-H instability problem is setup with a 2:1 density ratio and equal mass 



particles, identical to the setup described in [Price (2008) and using 512 x 512 particles in the low density 
fluid. Whilst shocks are smoothed by the artificial viscosity term (Figs. [8]and[9|, with only viscosity the 
jump in thermal energy at the contact discontinuity is not treated, resulting in a 'blip' in the pressure proflle 
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(Fig. [8|. This manifests in the 2D K-H problem as an 'artificial surface tension' effect, caused by the very 
same kind of pressure blip across the boundary (contact discontinuity) between the dense and the light 
fluids, suppressing mixing of the two. With conductivity applied (Fig. |9| the pressure is smooth across the 
contact discontinuity in both problems, which for the K-H problem means that the two fluids mix correctly. 
The lack of treatment of contact discontinuities in standard SPH codes (i.e., with no artificial conductivity 



term) explains the discrepancy between grid and SPH results somewhat infamously highlighted by Agertz 



et al.| ( |2007[ ). 

Fig. |10| shows the same shock tube in 2D (as already discussed briefiy in Fig. [5|. The main difference 
to ID is that the "noise" due to the particle resettling behind the shock front is visible (left subfigure). It 
is often asserted that SPH performs poorly on 2D shocks for this reason, however the noise can be very 
effectively minimised (at some additional cost) by employing the Mq quintic kernel instead of the cubic 
spline (right subfigure), giving results comparable to the ID version and illustrating in practice how the 
higher M„ kernels can be used to obtain convergence in SPH. 



7. Smoothed Particle Magnetohydrodynamics from a Lagrangian 

We can follow the same general approach to constructing an SPMHD algorithm as for hydrodynamics 
(Sec.|3]): Write down the Lagrangian, use appropriate physical constraints and use this to consistently derive 
the resultant equations of motion. 



7.1. MHD Lagrangian 



For MHD, the Lagrangian is given by (Price and Monaghan 2004b) 
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Ub{Pb, Sb) - 



2a*o Pb 



(106) 



corresponding simply to the subtraction of a magnetic energy term from the hydrodynamic verison (Eq. 16 ) 
In the continuum limit this corresponds to the standard MHD Lagrangian used by many authors (e.g 



Newcomb 1962 Field 1986) 
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(107) 



The difference to the hydrodynamic case is that, unlike the thermal energy term, neither the magnetic 
field B nor the change in the magnetic field can be written directly as a function of the particle coordinates, 
so we cannot straightforwardly employ the Euler-Lagrange equations (21). Instead, we can use the more 
general form of the variational principle given by 5S = J SLdt — (Sec. 3.2), where from (106) we have 
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(108) 



and the perturbation is with respect to a small change in the particle coordinates 5v. So we can derive the 
equations of motion provided that we are able to express the change in the magnetic field (5B as a function of 
the change in particle coordinates - equivalent to being able to write down an expression for the Lagrangian 
time derivative dB/dt [or equivalently, d(B/ p)/dt] since d/dt = 5/5t. In other words, in order to derive the 
SPMHD equations of motion it is necessary to specify not only the density estimate but also the manner in 
which the magnetic field is evolved. 



7.2. SPMHD formulation of the induction equation 

Given the induction equation for (ideal) MHD, written in the Lagrangian form 



d_ 
dt 



(109) 
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and using our general formulations of SPH derivatives given in Sec. 4.6 it is straightforward to write down 
an SPH version of the form 



dt \ Pa 



(110) 



which is equivalent to using the antisymmetric derivative (81 1) with <f> ^ p. 
1.3. Equations of motion 



The perturbations required in (108), from (56) and (110 1 are therefore given by 



Spb = ^^mdSn-Src) -^^bWbcihb) 
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(111) 
(112) 



giving, from (108) and using (25) 
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(113) 



where 5ba = Srb/Sva. Integrating the velocity term by parts, simplifying the double summations using the 
Kronecker deltas and the antisymmetry of the kernel gradient, and assuming the perturbations are 
arbitrary, we find that the SPMHD equations of motion are given by 
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In tensor notation these can be written more compactly in the form 
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(115) 



where S^^ is the MHD stress tensor, defined according to 

1 



S'^ = -\P 



2po 



B' 



5^3 + ^B'BK 
Mo 



(116) 



As for the hydrodynamic case (Sec. [s]) it is readily seen that the equations of motion conserve linear 
momentum exactly, due to the pairwise symmetry in the force. However, the MHD equations, unlike their 
hydrodynamic counterparts, do not exactly conserve angular momentum, since the perturbation to the 
magnetic field, (112) - and hence the anisotropic force term derived from it - is not invariant to rotations. 
It is interesting to note that the anisotropic magnetic force term in (114) derives entirely from the numerical 
representation of the induction equation (110), whilst the isotropic term derives purely from the magnetic 
energy term in the Lagrangian and the density perturbation (111 I. 
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7.-^. Energy equation 

The evolution equations for thermal energy and entropy - in the absence of dissipation - are identical 
to their hydrodynamic counterparts. The total energy evolution can be deduced from the Hamiltonian as 
in Sec. |3.4.2[ The corresponding expression for MHD is given by 
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Using (115), (35), (56) and (110 1, it can be shown that the total energy evolves according to 
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and is thus also conserved exactly. This further implies an evolution equation for the specific energy e 
u + ir—B'^/p of the form 
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7. 5. Interpretation of the Hamiltonian SPMHD equations 

Having used the SPH forms of the continuity and induction equations in the form 



d 
di 



dp 




'di ^ 


^ dx'' ' 






P ) 


p dxi ' 



(120) 
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it is straightforward to show (using 55) that our expressions (115) and (119) derived above are SPH repre- 
sentations of the MHD acceleration and energy equations in the form 
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where S"^^ is the MHD stress tensor given by (116). That we have explicitly used the induction equation to 



derive the equations of motion and energy is useful to our later discussion regarding what is a consistent 



formulation of monopole terms in the MHD equations (Sec. 10.1 1 



7.5.1. MHD Example 1: Advection of a current loop 

Our first MHD example demonstrates that some problems that prove very difficult for Eulerian schemes 
are almost trivial in a Lagrangian scheme such as SPMHD. The problem involves the advection of a loop 



of current across the computational domain, was introduced by Gardiner and Stone (2005) to test their 



Athena MHD code and presents a challenging problem for grid-based MHD schemes. The setup used here 



is identical to that in Rosswog and Price (2007) and we refer the reader to that paper (or the ndspmhd 
setup file) for full details. The current loop itself is given by the vector potential Az — Aq{R ~ x^ + y"^), 
with Aq set to give a weak field (plasma /3 of 2 x 10^) (here we compute and evolve the magnetic field, B). 
All particles in the domain (—1 < a; < 1, —0.5 < y < 0.5) are given a constant initial velocity along the box 
diagonal, with the magnitude set such that t = 1 represents one crossing of the domain. The results are 



shown in Fig. 11 at t = (left panel) and after 1000 (this is not a misprint!) crossings of the computational 
domain (right panel). In the absence of the explicit addition of resistivity terms, there is no change in the 
current or the magnetic energy and the advection is computed exactly. 



30 



Figure 11: Advection of a current loop across a 2D domain. Magnetic field lines are shown at t = (left panel) and after 1000 
crossings of the computational domain (right panel). In the absence of the explicit addition of resistivity terms, there is no 
change in the magnetic field and the advection is computed exactly in SPMHD, irrespective of the direction of propagation, 
numerical resolution or advection velocity. 



8. The tensile instability in MHD 



The rather large caveat to using the equations of motion derived in Sec. |7.3| is that in MHD the total 
stress can become negative, meaning that the force between particles can become attractive rather than 
repulsive and the equations will be unstable to the tensile instability discussed in Sec. |5.3| The regime of 
instability is evident if we consider the force in just one spatial dimension (and at constant h), given by 
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so the force will become attractive (along the field lines) w hen ^B'^/fj.p > P, i.e., w hen magnetic pressure 
exceeds gas pressure. A more detailed stability analysis (e.g. Phillips and Monaghan 1985 Morris 1996b|a 



B0rve et al. 2004 1 confirms that this is also the criterion for instability in more than one spatial dimension. 
Thus, whilst for particular applications that remain in the regime where magnetic pressure is smaller than 
gas pressure it is possible to use the conservative formulation (e.g. Dolag, Bartelmann, and Lesch 1999), 
in most cases stabilising the tensile instability is the first and most basic requirement for a stable SPMHD 
algorithm. 

The physical reason for the MHD instability is that the momentum-conserving force corresponds to 
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which differs from equations of motion written in terms of the Lorentz force (where J = VxB//io), 
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(126) 



by the "monopole term" in (125). This term gives a force directed along the magnetic field and proportional 
to the (numerically non-zero) divergence of the magnetic field, and is the source of the instability in the 
conservative formulation if it cannot be counteracted by pressure. As a result, Eq. |126| was used as the 
basis of a number of early SPMHD formulations in which the force was simply computed using standard 



curl operators such as ( 


54 


) or ( 


79 


) (e.g. Mcgficki 


1995 


Gouveia Dal Pino 


2001 


Hosking and Whitworth 


2004 


I. 



However, the poor conservation properties of such 
formulations means that MHD shocks are not well captured. 

Dealing with such 'source terms' is not an issue unique to SPMHD and requires careful consideration in 



all numerical MHD formulations (see Sec. 10.1 ). Of course, V • B = should be zero physically due to the 
non-existence of magnetic monopoles in the Universe, so both (126) and (125) are valid in the continuum 
limit. The problem is that the divergence term is not exactly zero numerically - so one is forced to make a 
choice. As noted by Toth (2000), it is also not a simple matter of enforcing the V - B = constraint in some 
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Figure 12: The tensile instability in SPMHD: The figures show the particle distribution (left panel in each figure) and x- 
component of the magnetic field (right panel in each figure) in the 2.5D circularly polarized Alfven wave test using the 
(unstable) conservative SPMHD force (left figure) and with a stable formulation (right figure), shown after 1 wave crossing 
time. Untreated, the MHD tensile instability results in a catastrophic clumping of particles along the magnetic field lines 
(left figure) which proceeds to destroy the calculation. Physically, it can be attributed to non-zero divergence terms in the 
conservative form of the MHD force, and can be stabilised by explicitly subtracting the unphysical 'source term' (right figure). 



way, since what is necessary to achieve a force that is both conservative and exactly perpendicular to B is 
to constrain V • B to be zero in the discretisation that the term appears in the force equatior^^ In SPMHD 
this is equivalent to requiring both exact derivatives and exact conservation which, as discussed in Sec. |4.7[ 
does not appear to be possible. 



.1. Fix 1: subtract a constant from the stress 



The original paper by Phillips and Monaghan ( 1985 ) proposed a simple fix involving a prior sweep over 



the particles to find the maximum (negative) stress, which would then simply be subtracted (as a constant) 
from the stress in the equations of motion, giving 
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(127) 



which conserves momentum but not total energy. The caveats are that there is a computational cost 
involved to compute 5^^^ and if this term is large it can lead to unphysical effects in the simulation. On 
the other hand, this is a simple technique that removes the instability and has relatively few side effects 
provided the correction is small. It is particularly useful if, for example, the simulation is dominated by 
large (constant) external stresses, whereby explicitly subtracting the external component of the stress (e.g. 
due to an externally imposed magnetic field) can serve to stabilise the formulation. 



3.2. Fix 2: use a more accurate but non- conservative gradient estimate in the anisotropic force 



Morris ( 1996b ) proposed a compromise approach whereby one retains conservative form in the isotropic 



part of the force, but adopts a more accurate but non-conservative derivative estimate (that vanishes when 



l°This was initially thought impossible by |T6th|2000l though [T6"th|2002| later showed such a discretisation could be achieved 
for grid codes. However, there does not appear to be an equivalent formulation in SPMHD, since the tensile instability occurs 
even in one dimension where the divergence constraint can be trivially enforced using = const, but dB^jdx 7^ in the 
force equation (c.f. Eq. |124|l. 
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the stress is constant) in the anisotropic term. Adapted for variable smoothing lengths, the resultant force 
is given by ( [Price and Monaghan 2005 Price 2010) 
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(128) 

where \/Wab = [VWab(/ia) + VM^ab(/ib)]/2. The 'Morris approach' does not exactly conserve momentum 
or total energy, but the errors due to non-conservation are in practice quite small even on strong shock 
tube problems (see Price[ 2004). Note, however, that retaining local conservation in the isotropic term is 
important for the reasons discussed in Sec. [5] and also in order to obtain the correct jump conditions on 
MHD shocks. The caveat to this approach is that the remaining non-conservation of momentum can become 
problematic, especially in simulations run for a long time period, where secular effects can accumulate and 
cause a loss of symmetry, and that it cannot be turned "off" in the otherwise stable weak-field regime. 



8.3. Fix 3: subtract the unphysical source term 



Tve et al. (2001) suggested an approach whereby conservation form is retained, but the monopole 



source term is explicitly subtracted. Extending their expression to incorporate variable smoothing length 
terms, the corrected equations of motion are given by 
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The corrective term with _B* = i?* is equivalent to the source term added to the momentum equation in 
the eight-wave Riemann solver (Powell et al. 1999) - with the key point from an SPH standpoint being 



that the discretisation used for the correction term is identical to that in which the divergence term occurs 



in the force equation. B0rve et al. (20041 also showed that it is not necessarily requir ed to always choose 
Bl = Bl in order to achieve stability, showing in particular that B^ ~ ^B^ is sufhcient. B0rve, Omang, and 
( 2006 ) give a more general method for setting B^ such that - amongst other things - no correction 



Trulsen 



is applied when the gas pressure exceeds the magnetic pressure (we refer the reader to their paper for more 
details). Whilst adding the correction term violates exact conservation of momentum and energy, it does so 



only insofar as the divergence term in (1291 is non-zero. Like the Morris approach, (129) is a good general 



method for removing the tensile instability in SPMHD - but has the advantage of being able to be switched 
ojf in the regime where the conservative formulation is stable. 



8.4. Other methods 

Finally, other methods have been proposed for dealing with the tensile instability in other physical 
problems, such as elastic dynamics. In particular [Monaghan] ([2OOOI suggested explicitly adding a short 



range repulsive force between particles to counteract the unphysical attractive force, similar to what would 
occur in atoms. [Price and Monaghan (2004a) initially applied this approach to the MHD equations and 



it was shown by Price (2004) to be equivalent to a modification of the kernel gradient in the anisotropic 
part of the MHD force (the second term in Eq. 114). However it proved problematic in highly compressible 
calculations where the smoothing lengths could vary dramatically (leading to a difficulty in defining "short 
range") and at large negative stress was found to result in large errors in the numerical sound speed (see 



B0rve et al. (20011 approaches 



Price[ 2004 1 . It was therefore abandoned by Price and Monaghan ( 2005 ) in favour of either the Morris or 



8.4.1. MHD Example 2: Circularly polarized Alfven wave in 2.5D 

An example of the tensile instability in practice is shown in Fig. |12[ which shows the particle distribution 
(left panels in each figure) and x-component of the magnetic field (right panels in each figure) in the 2.5D 
circularly polarized Alfven wave test from Toth (20001. The initial conditions (Toth, 2000 Price and 
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Monaghanl 120051) are p = 1, P = 0.1, vn ^ 0, Bn ^l,vj_ = Bj_= 0.1 sin (27rr|[) and = B 



= 0.1cos(27rr||) 
30° with respect 
1 / sin 9 with periodic boundary 



with 7 = 5/3 (where ry = x cos + y sin 9) , with the wave vector directed at an angle of 9 
to the a;- axis in the computational domain < x < l/cos0, < y < 
conditions. 

Using a conservative SPMHD formulation (left figure, shown at i = 1) whole lines of particles can be 
seen to have attracted each other in the direction of Py , which proceeds to destroy the calculation at later 
times. With either the Morris or B0rve et al. approach (right figure, also shown at i = 1) the particle 
arrangement - and hence the wave 



is stable and propagates correctly. 



9. Dissipation terms and shocks in SPMHD 

The second important issue for SPMHD is the formulation of dissipation terms appropriate for MHD 
shocks. We can follow the same general principle as in Sec. 6.3.2 starting with (102 1 as the basic equation, 
except that in MHD the energy variable is given by 
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and the signal velocity in the kinetic energy term is generalised to 
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where v is the fastest wave speed for MHD (i.e., the fast MHD mode), given by 
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The key constraint in deriving dissipative terms in the other equations is that the contribution to the 
entropy from the dissipative terms in the total energy equation must be positive definite in order to satisfy 
the second law of thermodynamics. For the kinetic energy term this is satisfied provided a term is added 
to the acceleration equation (i.e., Eq. 101 1, giving a positive definite contribution to the thermal energy 
(Eq.[l04|. It may also be shown that the thermal energy term results in a positive definite contribution to 
the entropy (Appendix B of Price and Monaghan||2004a I . For the magnetic energy term, positivity means 
that the overall dissipative contribution to the thermal energy equation is given by 
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which is positive definite because Fah is negative definite (Fig. [2| and the terms inside the brackets are 
positive. Using (133), we can deduce the term which must arise in the induction equation using 
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giving a term in the induction equation of the form 
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Figure 13: A "1.75D" MHD shock tube problem, showing the formation of 7 discontinuities: A slow, Alfvcn and fast wave 
propagating in each direction, plus the contact discontinuity. Capture of all of these discontinuities requires the application of 
artificial viscosity, thermal conductivity and resistivity to treat jumps in v, u and B respectively. The M4 cubic spline kernel 
has been used. 



The interpretation of this term is straightforward using ( 95 ) : We have derived an artificial resistivity term 



dt 
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(137) 



For the artificial resistivity it is also better to use a signal velocity that differs from the one used in the 
viscosity. A simple choice (used in ndspmhd) is to use an averaged Alfven speed, e.g. 



1 



"A,b' 



(138) 



Price and Monaghan ( 2004a 2005 1 also show how an artificial resistivity can be derived that uses only 



components of the magnetic field perpendicular to the line of sight. However, Price and Monaghan (20051 



found that the version using the jump in total energy (Eqs. 130 135 ) performed better in practice. The 
parameter as can also be evolved using a switch similar to ( |105[ ), where Price and Monaghan ( 2005 ) suggest 
a source term of the form 

5 = max(^,M,. (139) 



Note that artificial thermal conductivity and resistivity should in general be applied regardless of whether 
or not the particle pair is approaching or receding. 

9.1. MHD Example 3: shock tube problems in MHD 

The application of specific dissipative terms for shock-capturing is illustrated in Fig. [13) which shows a 
"1.75D" MHD shock tube problem (i.e., ID but with 3D magnetic and velocity fields) from Ryu and Jones 
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Figure 14: The Brio-Wu MHD shock tube problem in 2D, computed using 800 X 24 particles and the Mg quintic kernel. With 
the quintic, results comparable to the ID solution can be obtained. The numerical solution from [Balsara] | |1998| is given by 
the solid (red) line. All particles are plotted. 



(1995), using 800 particles and the standard M4 cubic spline kernel. The setup is (p, P, u 



X ; 1 ^2 ; Py ; 



[r0870.95, 1.2, 0.01, 0.5, 3.6/(47r)i/2,2/(47r)i/2] f^^. ^ ^^^^^^ for x > we have (p, P,v^,Vy,v^, By, B^) = 
[l,l,0,0,0,4/(47r)i/2,2/(47r)i/2j^ with B^ = 2/(47r)i/2 and 7 = 5/3. The main point is that this shows 
all 7 possible discontinuities in MHD: Three waves (slow, Alfven and fast modes) propagating in each 
direction, plus the contact discontinuity. Since each requires treatment, we require artificial viscosity, thermal 
conductivity and resistivity to treat the jumps in velocity, thermal energy and magnetic field respectively. 
Fig. 



14 shows a 2D version of the 
for MHD codes (e.g. [Stone et al.] |1992 
using 800 X 24 particles and the Me 
{p,P,v^,Vy,By) = [0.125,0.1,0,0,- 
SPMHD particles are plotted at t 



Brio and Wu (1988) shock tube test, considered a standard test 



Dai and Woodward 1994 Ryu and Jones 1995 Balsara 1998) 



quintic kernel. Initially {p,P,Vx, 



'Jy, By) 



[1,1,0,0,1] for a; < 0, with 



1] for a; > with B^ = 0.75 everywhere and 7 = 2.0. Profiles of all 
0.1 whilst the numerical solution from Balsara] ( 1998 ) is given by the 
solid (red) line. The solution is comparable to that achievable in ID (e.g. Rosswog and Price [2007 ) - at 
this resolution the only noticeable deviation from the Balsara solution is the slight mismatch in thermal 
energy at x « 0.15-0.3, a result of the low spatial resolution (lOx lower than in the a; < region) in this 
low density region. 



10. The divergence constraint in SPMHD 

The third major issue, common to all numerical MHD codes, is how to enforce the V • B = constraint. 
The problem arises because the constraint occurs only as an initial condition in the MHD equations, since 
from the induction equation we have 

— =-Vx(vxB); _(V.B) = 0, (140) 

The problem is that the truncation error in the numerical solution means that such a condition cannot 
be maintained indefinitely. What to do about it falls into three general approaches: "ignore" , "clean" or 
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"prevent" the numerical growth of monopole terms. We will discuss each of these in the context of SPMHD, 
but first turn to the issue of formulating consistent equations given that V • B 7^ 0. 

10.1. V • B 7^ 0; Monopole terms in the evolution equations 

We have already touched on the 'source terms' issue in the context of the tensile instability in SPMHD, 
which arises from the fact that the MHD force written as the divergence of a stress tensor differs from the 
exactly-perpendicular-to-B Lorentz force by a term proportional to V • B - and that this term does not 
vanish even in the trivial case of ID where = const. A similar issue arises in the induction equation 
which in "conservative form" is given by 



(140) 
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l)t 



= V • (vB - Bv) = -(v • V)B + (B • V)v - B(V • v) + v(V • B), 



(141) 



identical to (1401, and conserving the volume integral of the flux J BdV. However, if the monopole term 



v(V • B) is neglected, then we have (using the Lagrangian time derivative) 

dB 



= (B- V)v-B(V-v) 



(142) 



taking the divergence of which shows that the monopoles evolve according to an equation similar to the 
continuity equation for density 

|^(V-B) + V-(vV-B), (143) 
at 

meaning that the volume integral f^^(V • B)dV, and by Green's theorem the surface integral of the flux, 
§g B • dS is conserved. So adopting ( 142 ) represents a 'monopole conserving' form of the induction equation. 



such that the surface flux integral will be conserved even if the divergence errors are non-zero. On the basis 



of this idea Powell et al. (19991 suggested the 8- wave formulation, with the induction equation given by 



( 142 ) and the momentum and energy equations given by 
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However, Janhunen (2000) and Dellar (2001) have strongly argued that a consistent formulation in the 



presence of monopoles, whilst adopting (142 1 instead of (140), should nevertheless still conserve linear (but 
not angular) momentum and energy. For SPMHD we have already proved that this is the case because 
we derived the momentum-conserving force (115) from the surface-flux-and- monopole conserving induction 
equation ( |109[ ). This is a somewhat moot point in practice though, since we require the subtraction of the 

giving 



.3) 



source term in the momentum equation anyway in order to stabilise the tensile instability (Sec. 
essentially the 8-wave formulation. 

10.2. V • B w 0; The "ignore" approach 

The simplest approach to the divergence constraint is to do nothing at all. That is, simply monitor 
the divergence error and "trust" the simulations if it remains small. In SPMHD it is usual to monitor the 
divergence error using the dimensionless quantity 



W - B 



(146) 



typically hoping that it remains smaller than a few percent. For a number of problems, including many of 
the test problems shown in this paper, this approach works surprisingly well. There is also the question 
of what one means by the "divergence error" , since ( 146 ) can be computed using a variety of operators 



(Sec. 4.4) and the 'stencil' in which one measures divergence changes in SPH as soon as one changes the 
particle positions. However there are clearly many problems where this approach is not sufficient. 
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10.3. V • B — !■ 0; The "clean" approach 

Divergence cleaning in SPMHD also faces the problem of defining what V-B = actually means. Clearly 
from a stability point of view it would be desirable for the divergence to be zero in the discretisation in 
which it enters the force equation, but we have already seen that this is not possible even in ID because 
the conservative SPH derivative does not vanish for constant functions. However, it may be hoped that 
cleaning in a particular discretisation will at least produce solenoidal fields to the order of the truncation 
error [i.e., ©(/i^)]. Approaches to divergence cleaning fall into three general categories based on the solution 
to an elliptic, parabolic or hyperbolic partial differential equation. 

10.3.1. Elliptic cleaning 

Elliptic cleaning corresponds to using projection methods: Solving either for a correction term, 



B = B* - V(/); 



= V B*, 



or alternatively for the physical (solenoidal) component of the field, 

B = VxA; V^A = -VxB* 



(147) 



(148) 



SPH methods for both of the above have been discussed by Price and Monaghan (2005). The main issue is 



that, if one simply uses the Green's function solution - analogous to the computation of gravitational forces 
- with the source function computed with a chosen div or curl operator, the projection is only approximate 
(that is, the Poisson equation is not discretely satisfiec|^. For 3D simulations with particles on individual 
timesteps, it is also impractical and highly inefficient to compute a global projection every n timesteps. 

10.3.2. Parabolic cleaning 

A parabolic approach corresponds to adding a V • B diffusion term to the induction equation, i.e.. 



~dtj 



(149) 



clean 



Having formulated our artificial dissipative terms using the jump in total magnetic energy ( 130 1 means that 



the artificial resistivity already contains a term of this form (see 1371, so applying artificial resistivity (with 



a switch that responds to the divergence, Eq. 139 ) corresponds to a form of parabolic divergence cleaning. 



This is partly why the "ignore" approach is not as bad as it sounds, and has successfully used for a number 



of problems (e.g. Dolag and Stasyszyn 2009 Kotarba et al. 2010 Biirzle et al. 2010 ) 



10.3.3. Hyperbolic/parabolic cleaning 



Dedner et al. (|2002|) suggested a hyperbolic cleaning method for the MHD equations, given by 

-B(V-v) + (B- V)v- VV-, 



dB 

lit 



where -0 is an additional variable that evolves according to ( Price and Monaghan 2005 ) 

-4(V-B) 



dip 
~dt 



(150) 



(151) 



where for SPH it is natural to use Ch 



:ig and T = ahh/vsig where an € [0, 1] is a dimensionless parameter. 
The first term in ( 151 1 corresponds to hyperbolic propagation of the "divergence wave" (at speed Ch), whilst 



^^The approximate nature of the projection is related to the question of how one should "soften" the Green's function 
solution, similar t o the case for gravity. What is required is a generalisation of the softening formulation given by [Price and] 
|Monaghan| | |2007[ l which shows how the gravitational Poisson equation can be discretely satisfied if the SPH density is used as 
the source term, so this would be an approach worth pursuing. 
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Vj/t is a parabolic decay term (with decay time r). This can be shown combining the equations (assuming 
V = const) to yield a wave equation, 



1 d^^p 



d4> 
'~dt 



= 0. 



(152) 



As reported by Price and Monaghan (20051 and others (e.g. Mignone and Tzeferacos 2010), the original 



Dedner et al. (20021 formulation contained a parameter that is not dimensionless, but that is equivalent 
to the wavelength of critical damping Ac in the parabolic decay term (evident by writing r = Xc/ch in 



the above). Indeed, Price and Monaghanj ( 2005| found that whilst using ( 150 l-( |l5T ) could be effective at 
cleaning well-resolved divergence errors (Ac ^ h), it was found to have little or no effect when used in "real" 
problems (where Ac ^ ft.), giving at most a factor of ^ 2 reduction in V • Bmax- It was also found that a 
poor choice of parameters could increase the divergence error substantially via constructive interference of 
divergence waves. 

More recently Stasyszyn & Dolag (following Dolag and Stasyszyn 2009 1 report better success using this 
method, by two adaptations: i) accounting for the magnetic energy dissipation in the energy equation, and 
ii) implementation of a limiter to control the amount by which the magnetic field can change in a single 
timestep (Dolag 2010, priv. commun.). Nevertheless, it does not appear that the method is being used on 
real problems (e.g. Kotarba et al. 20101, suggesting that it remains ineffective. 



10.4- V • B = 0; The "prevent" approach 

The third approach is to try to 'prevent' divergence errors altogether by formulating the MHD equations 
such that V • B = is satisfied by construction. Unfortunately, the approach most successfully adopted 
in grid based schemes - the constrained transport method (Evans and Hawley 1988 Gardiner and Stone 



2005 ) - cannot be directly applied in an SPH context because it requires the computation of surface rather 



than volume integrals. 



10.4.1. Euler potentials 

Perhaps the closest to a 'constrained transport' approach in SPMHD is the formulation in terms of Euler 
potentials (Stern 1970 1976) (also referred to as "Clebsch variables", e.g. by Phillips and Monaghan|1985 ), 
aE and Pe, where the magnetic field is expressed as 



B = Vue X 



(153) 



Taking the divergence shows that V • B is satisfied by construction. A further advantage is that for 



ideal MHD, the induction equation (109) becomes simply 



daE 
dt 



-0, 



df3E 
dt 



= 0, 



(154) 



1966) 



which corresponds physically to the advection of magnetic field lines by Lagrangian particles ( Stern 
That is, the Euler potentials refiect the physical fact that the field lines are 'frozen' to the fiuid in ideal 
MHD and are therefore a natural description for a Lagrangian scheme. For SPH it is straightforward to write 



down expressions for (153) using any of our standard first derivative operators (Sec. 4.4) and to simply use 
the resultant B in the SPMHD equations of motion - though this will give only approximate conservation 
of energy, see Price |2010. However, in order to capture shocks, some dissipation in the magnetic field is 



required. For this reason Price and Bate ( 2007 1 and Rosswog and Price ( 2007 ) introduced dissipation terms 
of the form 



da% 
dt 



diss 



E 



rub B 
—asv 

Pab 



sig 



) Fab\ 



dP% 

dt 



diss 



E 



TTlb B 

— asv 

Pab 



sig 



{P% - P'k) Fab, (155) 



corresponding to the continuum equations (assuming r/M ~ const) 

dpE 



daE „2 
^ - ^mV a- 



dt 



(156) 
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where rjM is given by (137). However it is clear that this is not a consistent formulation of (physical) 



resistivity for the Euler potentials, since we have neglected mixed terms that are of comparable magnitude 
to those in (156) (e.g . Brandenburg 2010). It is also difhcult to ensure that the contribution to the entropy 
is positive definitern although Price (20101 shows how the latter requirement can be achieved for the vector 



potential, and similar formulations should be used in place of (155) for the Euler potentials 



Nevertheless, the control of the divergence errors means that the Euler potentials have found successful 



application to a range of problems (e.g. Price and Bate 2007 Dobbs and Price 2008 Kotarba et al. 2009 ) 



most notably in star formation (e.g. Price and Bate 2007 2008 , 2009 ) since one is able to form stars in the 
presence of magnetic fields without blow-up of the numerical solution due to divergence errorj^ The caveat 
is that the Euler potentials do not capture the full physics of MHD, since the Helicity A • B is identically 
zero, meaning that topologically non-trivial fields cannot be represented (or by corollary, created during a 



simulation) using (153)-(154). For example, the Euler potentials cannot be used to follow multiple complete 



windings of a magnetic field in a rotating fluid, a corollary of the fact that the potentials are simply advected 



by the particles (Eq. 154), so reconstruction of the field via (153) requires a one-to-one mapping between 
the initial and final particle positions. These restrictions are easily demonstrated by simple test problems 



( Brandenburg 2010 ) and mean in practice that one is limited to studying problems with initially simple field 



geometries, and that important field winding processes (i.e., any dynamo action) are missed. On the other 
hand, it may be possible to formulate a more general Euler potentials description without such limitations. 

10.4-2. The vector potential 

Another possibility for a divergence-free formulation without the restrictions of the Euler potentials is 
to simply employ the vector potential B = V x A directly. Indeed in 2D they are exactly equivalent, since 
for a 2D field = Az{x, y) and = z. The main disadvantage of the vector potential in 3D is that the 
evolution equation is complicated and requires an appropriate gauge choice. For example, setting the scalar 
potential = v • A one can obtain a Galilean-invariant evolution equation for A of the form (Price 20101 

dA 



dt 



= A X (V X v) - (A • V)v + V X Bej^f = -A^Vv^ 



X Bp 



(157) 



(106). In particular, Price (20101 shows that by writing the curl according to 



Using this, one can proceed to derive an SPMHD algorithm for the vector potential from the Lagrangian 

(158) 



B<j = (V X A)a + Be^t = — "^6(Aa - Af,) X WaWab + B^^t, 

Pa V 



the effective evolution equation for B is given by 

dBa 

Pa 

( dAa dAb 



dt 



+ 



- ^ mb (A„ - Afc) X [(v„ - Vb) • V] VaWa, 



V dt 



V7 mr Bint dpa 

X VaWab -rr- 

dt I Pa dt 



Thus, writing the evolution equation for A, (157), in the form 

dA A^ 

— ^ = ^ X! "^^^^i ~ '^D'^oWah + (V X Be^t)a, 
Pa , 



(159) 



(160) 



^ ^Rosswog and Price | |2007| simply used Eq |l33| with B calculated from l |l53| such that the contribution is guaranteed to 
be positive. However, t he total energy is only approximately conserved by this procedure. See|Price| (j2010). 

^ ^Price and Rosswog| | [2"006| also employed the Euler potentials to simulate magnetic field growth during JNeutron star mergers. 
Whilst similar results were found using a B-based formulation, the simulations with Euler potentials were later found to be 
incorrect, showing exaggerated field growth because the boundary conditions on the potentials for the stars were not correctly 
accounted for. 
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one can self-consistently derive equations of motion for the vector potential from ( 108 ) by writing the 
perturbation to 6'B in terms of 6 A using (159) and 5 A in terms of Sr using (160), giving 



dt 




3_ n2 \ 

''''' V,M^,, 



Pi 



[(A, - Afc) X V] \ V,PF„, 



2 v-^ /Bq Bfc 



Mo 



b 



f, V Pa Pb 



MO V 



B„ 



B. 



Pa Pb 



Bea;t ' aWab 



Pi 



Pb 



where J the magnetic current is computed using the symmetric curl operator (the conjugate to 1581 

T _ (V X B)„ 



Mo 



Ma 

Mo 



E 



rrib 



Ba 

Pi 



Bb 

Pi 



(161) 



(162) 



The proof that the above equations do indeed translate into an continuum expression of the conservative 
MHD equations, as well as a more general version taking full account of smoothing length gradient terms is 



given by Price (20101. However, the major flaw, from simple inspection of (161), is that the isotropic (first) 
term will give a negative pressure - and thus the tensile instability - whenever 3-B^/(2/io) > P, a much 
more restrictive regime than for the usual SPMHD equations (see Sec.js]). This is verified in test problems, 
but unlike the standard SPMHD instability, proves difficult to stabilise without significantly affecting the 
solution. A second instability was also found to arise purely due to the evolution of A according to (157), 



resulting in unconstrained growth of vector potential components, most likely related to the need to explicitly 
enforce a gauge condition alongside the evolution of the vector potential. Price (2010) therefore concluded 
that using the vector potential was not a viable approach for SPMHD. 

11. Summary 

In summary, we have given an overview of SPH methodology, starting with the density estimate as the 
basis of all SPH formulations (Sec. [2]) and showing how the equations of motion and energy can be self- 
consistently derived from the density estimate using a variational principle (Sec. [3|. Kernel interpolation 
theory has been introduced mainly as a way of interpreting the SPH equations, and we have discussed how 
linear error analysis can be used to construct more accurate and very general derivative operators (Sec.|4]). 
In Sec. [5] the importance of local conservation was highlighted with respect to maintaining a regular particle 
distribution and thus accurate derivatives, giving us an understanding of how the tensile instability arises in 
SPMHD and why one should be careful in setting the ratio of smoothing length to particle spacing. Second 
derivatives in SPH were discussed in Sec. |6j mainly as a way of formulating dissipative terms necessary 
to capture shocks and other kinds of discontinuities. In the second half of the paper, we have shown how 
SPMHD, like SPH, can also be formulated from a variational principle (Sec. [7| and have addressed the 
three main issues with regards to the accuracy of SPMHD: The tensile instability (Sec.|8|, the formulation 
of shock-capturing dissipation terms (Sec. |9| and the enforcement of the divergence-free condition on the 
magnetic field (Sec. 10). Finally, this paper marks the public release of the NDSPMHD SPH/SPMHD code 
that can be used to test and verify all of the ideas and methods that have been discussed. 
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